error in qpsk modulation

Hi....

i am trying to implement a mimo-ofdm systems in which i have selected qpsk as a modulation technique. at the transmitter i am getting a proper constellation diagram. i have implemented the receiver section in which i am not able to properly get the constellation diagram for the same. please find attached code .

clc;
clear all;
close all;

% All Parameters Initialization  to be done here - begin
M=4;%  4,8,16,32,64,128,256,512,1024- 2,3,4,5,6,7,8,9,10 size of constellation 
k=log2(M); % number of bits/symbolexit
number_of_data_points=64;         % FFT Size or total number of subcarriers=N=64= IFFT size
numPackets=number_of_data_points;
block_size=16;                    % Block Size, no of rows,data subcarriers
BW=8*10^6;                       % Bandwidth of the OFDM system
Nsp=52;
%Derived Parameters
deltaf=BW/number_of_data_points;%bandwidth of each subcarriers including used and unused subacarriers
TFFT=1/deltaf;                  % ofdm_symbol_duration
TGI=TFFT/4;                     % duration of the cyclic prefix= guard interval duration 25% of the subcarriers to be used as a GI
% possible guard band values 1/4, 1/8, 1/16 and 1/32
Tsignal= TGI+TFFT ;             % total ofdm symbol duration
Ncp=number_of_data_points*(TGI/TFFT); % length of the cyclic prefix
cp_len=ceil(0.25*block_size);
fft_points=number_of_data_points;
ifft_points=number_of_data_points;
Nt=2;
Nr=2;
cp_start=block_size-cp_len;         % with block_size=4,1x1;  block_size=2,1
cp_end=block_size;                  % with block_size=4,1x1;with block_size=2,1x1,2
actual_cp=zeros(cp_len, block_size);%1x4; with block_size =2, actual_cp= 1x2;
data_subcarriers=cp_start+Nsp;
nsf=8/7;
Fs=ceil(nsf*BW);
EbNo=0:10;
% initialization of SNR 
SNRstart=1;
SNRincrement=1;
SNRend=9;
index=0;
EbN0=0:10;   % Eb/N0 
E12=1; 
SNR=10^(EbN0(k)/10);                     %  Eb/N0(dB) to linear value
standar_dev =sqrt((E12/(2*SNR))*(15/11));  % Standard Deviation of Noise
bit_error=0; 


%-------------------------------------------------------------------------%
%1.  1st Block  of user 1 -Implementation of binary data                             %
%1st Block -Implementation of bianry data- generation                     %
%-------------------------------------------------------------------------%
x1U=randi([0 1],1,64);

figure (1)  %data input plotbr=BW; %Let us transmission bit rate  1000000
stairs(x1U(1,:), 'linewidth', 3)
xlabel('Data Points')
ylabel('Amplitude')
title('Information before Transmiting');
axis([ 0 63 -1 2 ]);
%-------------------------------------------------------------------------%
c=64;%size(x1,2);
n=64;%size(x1,2);% number of information bits/ message bits/ length of message vector...
wc=4;         % enter the number of ones in each columns
wr=(n.*wc)./c; %number of ones in each row
H1=zeros(c,n); %generate the empty matrix and start assigining the 1's c=n; n=k
j=1;
jj=wr;

for i=1:c./wc
H1(i,j:jj)  =1;
j=j+wr;
jj=(i.*wr)+wr;
end
for i=1:wc-1
    for ii=1:n
        colind=(round(rand(1)*(n-1)))+1;
        rCol=H1(1:c./wc,colind);
        H1((i.*(c./wc))+1:((i.*(c./wc))+1+(c./wc)-1),ii)=rCol;
    end
end

br=BW; %Let us transmission bit rate  1000000
fc=br; % minimum carrier frequency
Tc=1/br; % bit duration
tmod=Tc/99:Tc/99:Tc; % Time vector for one bit information

%-------------------------------------------------------------------------%
%------- encoded data by using LDPC Encoder and Parity Check Matrix-------%
%-------------------------------------------------------------------------%

data_encoded1=x1U*H1;

mod_data1=mod(data_encoded1,2);
figure (2)   % data encoded by using LDPC 
stairs(mod_data1(1,:), 'linewidth', 3)
% stem(mod_data, 'linewidth',3), grid on;
title('  Information Encoded by LDPC Encoder ');
xlabel('Data Points')
ylabel('Amplitude')
axis([ 0 63 -1 1.5]);


%-------------------------------------------------------------------------%
%---------3. modulation of the data by using psk------------------------- %
%-------------------------------------------------------------------------%

data_NRZ1=2*mod_data1-1; % Data Represented at NZR form for QPSK modulation
s_p_data1=reshape(data_NRZ1,2,length(x1U)/2);  % S/P convertion of data
ymod1=[];
ycomplex1=[];
y_in1=[];
y_qd1=[];
% for jmodsp=1:size(s_p_data,1)
    for imod = 1:length(x1U)/2
    y11 = s_p_data1(1,imod)*cos(2*pi*fc*tmod); % inphase component
    y21 = s_p_data1(2,imod)*sin(2*pi*fc*tmod);% Quadrature component
    y_in1 = [y_in1 y11]; % inphase signal vector
%     y_in1=reshape(y_in,16,[]);
    y_qd1 = [y_qd1 y21]; %quadrature signal vector
    ymod1 = [ymod1 y11+y21]; % modulated signal vector
    ycomplex1 = [ycomplex1 (sign(y11))+(1i*sign(y21))];
    end
% end

psk_mod1=ycomplex1;
% psk_mod=ymod;
Tx_sig11=ymod1; % transmitting signal after modulation
Tx_sig1=ycomplex1; % transmitting signal after modulation
Tx_sig_complex1=ycomplex1;
% figure (3)
scatterplot(Tx_sig_complex1)
ttmod=Tc/99:Tc/99:(Tc*length(x1U))/2;

figure (4)
subplot(2,1,1) % plot of inphase compnent
stairs(s_p_data1(1,:),'linewidth',1.5),grid on;
axis([ 0 31 -2 2])
subplot(2,1,2)
plot(ttmod,y_in1,'linewidth',1.5), grid on;
title(' wave form for inphase component in QPSK modulation ');
xlabel('time(sec)');
ylabel(' amplitude');

figure (5)
subplot(2,1,1) % plot of inphase compnent
stairs(s_p_data1(2,:),'linewidth',1.5),grid on;
axis([ 0 31 -2 2])
subplot(2,1,2)
plot(ttmod,y_qd1,'linewidth',1.5), grid on;
title(' wave form for Quadrature component in QPSK modulation ');
xlabel('time(sec)');
ylabel(' amplitude');

figure (6)
subplot(3,1,1);
plot(ttmod,y_in1,'linewidth',1.5), grid on;
title(' Wave form for In-phase component in QPSK modulation ');
xlabel('Time(sec)');
ylabel(' Amplitude');

subplot(3,1,2);
plot(ttmod,y_qd1,'linewidth',1.5), grid on;
title(' Wave form for Quadrature component in QPSK modulation ');
xlabel('Time(sec)');
ylabel(' Amplitude');

subplot(3,1,3);
plot(ttmod,Tx_sig1,'r','linewidth',1.5), grid on;
title('QPSK modulated signal (sum of inphase and Quadrature phase signal)');
xlabel('Time(sec)');
ylabel(' Amplitude');
% to  find out the number of columns after reshaping
num_colums1=ceil((size(Tx_sig_complex1,1)*size(Tx_sig_complex1,2))/block_size);

%-------------------------------------------------------------------------%
%              5.to perform serial to parallel conversion                 %
%-------------------------------------------------------------------------%
% reshaped_mod_data = reshape(Tx_sig_complex,[block_size,num_colums]);
reshaped_mod_data1 = reshape(Tx_sig11,[block_size, num_colums1]);
% data divided into 4 subcarriers
reshaped_mod_data_11=reshaped_mod_data1(:,1);
reshaped_mod_data_21=reshaped_mod_data1(:,2);
reshaped_mod_data_31=reshaped_mod_data1(:,3);
reshaped_mod_data_41=reshaped_mod_data1(:,4);

%-------------------------------------------------------------------------%
%            6.     To Perform IFFT and Cyclic Prefix                     %
%                 to create empty matrix to put ifft data                 %
% 7.operate column wise and do cyclic prefix, column wise cyclic prefix is done
%---------------------------+----------------------------------------------%

% ifft_data_matrix=zeros(size(reshaped_mod_data,1),size(reshaped_mod_data,2));

for ireshape=1:num_colums1
    ifft_data_matrix1(:,ireshape)=ifft(reshaped_mod_data1(:,ireshape)); 
    for jcplen= 1:cp_len  % compute cyclic prefix data
        actual_cp1(jcplen,ireshape) = ifft_data_matrix1(jcplen+cp_start,ireshape).';
    end
    ifft_data_after_cp1(:,ireshape) = vertcat(actual_cp1(:,ireshape),ifft_data_matrix1(:,ireshape));% perform cyclic prefix
end

%-------------------------------------------------------------------------%
%           8.    to perform parallel to serial conversion                %            
%-------------------------------------------------------------------------%
parallel_serial_tx1 = reshape(ifft_data_after_cp1(:,:),size(ifft_data_after_cp1(:,:),2),size(ifft_data_after_cp1(:,:),1));


[numrows,numcols]=size(parallel_serial_tx1);

x=parallel_serial_tx1;

snr=0:2:20;

error1 = zeros(1, numPackets); BER1 = zeros(1, length(snr));


Taps=1;                                         % Number of Taps 
p1=0.5/2.3;                                     % Power of Tap1 

gain1=sqrt(p1/2)*[randn(1,numrows) + j*randn(1,numrows)];   % Gain for Tap1 


Nc=16;
cp=4;
x11=x(:); 

x12=reshape(x11,1,length(x11)); 
i=1:length(x12);         
delay1=1;  
for i=delay1+1:length(x12) % Producing one sample delay in Tap2 w.r.t. Tap1 
   x13(i)=x(i-delay1); 
end 

x1=reshape(x13,(Nc+cp),ceil(length(x13)/(Nc+cp))); 
ch1=repmat(gain1,numcols,1)';      


data_channel=x.*ch1;  % Passing data through channel  

%------------------------Addition of AWGN noise --------------------------- 
data_noise1=data_channel(:); 
data_noise2=reshape(data_noise1,1,length(data_noise1)); 
noise = 1/sqrt(2)*[randn(1,length(data_noise2)) + j*randn(1,length(data_noise2))];

snr=0:2:20;
for i = 1:length(snr) 
    y = data_noise2 + (sqrt(1)*10^(-snr(i)/20))*noise; %Addition of Noise  
end 
%--------------------------Receiver --------------------------------------- 
 
data_received =y;           %fadded data received with awgn noise 
 

% serial_parallel_form_rx = reshape(received_signal,[ifft_rows,ifft_colums]).';%64x2
serial_parallel_form_rx = reshape(data_received,[numrows,numcols]).';%64x2

% To Remove Cyclic Prefix
serial_parallel_form_rx(1:cp_len,:,:)=[];

% To Perform FFT operation
receivd_signal_fft_rx=fft(serial_parallel_form_rx);

%  To Perform Parallel to Serial Conversion
 parallel_serial_form_rx= reshape(receivd_signal_fft_rx,1,[]);
 
 
%  To pass through a Demodulator
Rx_data=[];
Rx_sig=parallel_serial_form_rx; % Received signal
Rx_sig_real=(real(Rx_sig));
Rx_sig_imag=(imag(Rx_sig));


for (idemod=1:1:length(x1U)/2)
    %%XXXXXX inphase coherent dector XXXXXXX
     Z_in=Rx_sig_real((idemod-1)*length(tmod)+1:idemod*length(tmod)).*cos(2*pi*fc*tmod);   
    % above line indicat multiplication of received & inphase carred signal
    
    Z_in_intg=(trapz(tmod,Z_in))*(2/Tc);% integration using trapizoidal rule
    if(Z_in_intg>0) % Decession Maker
       Rx_in_data=1;
    else
       Rx_in_data=0; 
    end
    
    %%XXXXXX Quadrature coherent dector XXXXXX
    Z_qd=Rx_sig_imag(((idemod-1)*length(tmod))+1:idemod*length(tmod)).*sin(2*pi*fc*tmod);
    %above line indicat multiplication ofreceived & Quadphase carred signal
    
    Z_qd_intg=(trapz(tmod,Z_qd))*(2/Tc);%integration using trapizoidal rule
    if (Z_qd_intg>0)% Decession Maker
        Rx_qd_data=1;
    else
        Rx_qd_data=0;
    end
    Rx_data=[Rx_data  Rx_in_data  Rx_qd_data]; % Received Data vector
end
 
psk_demod_mod=Rx_data;
Hinv=H1.';
[N1decode,N2decode]=size(H1);
iter=3;
for idecode=1:iter
    %     for jpr=1:Nt
    for jd=1:N1decode %  for loop along row 
        cid=find(H1(jd,:));
        for kdec=1:length(cid)
            E(jd,cid(kdec))=mod(sum(psk_demod_mod(cid))+psk_demod_mod(cid(kdec)),2);
        end
    end
% end
%     for jpc=1:Nt
for jd=1:N2decode %for loop  along columns
    rid=find(H1(:,jd));
    numberofones=length(find(E(rid,jd)));
    numberofzeros=length(rid)-numberofones;
    if(numberofones==numberofzeros)
        ydd(jd)=psk_demod_mod(jd);
    elseif(numberofones>numberofzeros)
        ydd(jd)=1;
    elseif(numberofones<numberofzeros)
        ydd(jd)=0;
    end
end
 ydec=ydd;
end
rece_decoder=psk_demod_mod*Hinv;
rece_decoder_mod=mod(rece_decoder,2);


figure (7)
subplot(3,1,1);
stairs(Z_in,'linewidth',1.5), grid on;
title(' Wave form for In-phase component in QPSK modulation ');
xlabel('Time(sec)');
ylabel(' Amplitude');
subplot(3,1,2);
stairs(Z_qd,'linewidth',1.5), grid on;
title(' Wave form for Quadrature component in QPSK modulation ');
xlabel('Time(sec)');
ylabel(' Amplitude');
subplot(3,1,3);
stairs(Rx_data,'r','linewidth',1.5), grid on;
title('QPSK modulated signal (sum of inphase and Quadrature phase signal)');
xlabel('Time(sec)');
ylabel(' Amplitude');