vblast_qpsk.m

毕业设计matlab程序 VBLAST的ZF,MMSE,QR,ZF_IC,MMSE+IC译码程序

close all;
clear
echo off
Tx_n = 4;
Rx_n = 4;
index = 2;
%frame_length=2000000;
frame_length=1200;
SNRindB=0:2:20;

%***********************************************
%correlation
Rt=eye(Tx_n);
Rr=eye(Rx_n);
        for k=1:Tx_n
            for l=1:Tx_n
               if k>l
                  Rt(l,k)=0;
                  Rt(k,l)=Rt(l,k);
               end
            end
        end
        for x=1:Rx_n
            for y=1:Rx_n
              if x>y
                  Rr(x,y)=0;
                  Rr(y,x)=Rr(x,y);
               end
            end
        end             %不如cholesky分解效率高，但是可以随着天线数的变换方便改
%************************************************

for i=1:length (SNRindB),
    SNR(i)=10^(SNRindB(i)/10);
    a=rand(1,frame_length);
    for L=1:frame_length
        if a(L)>.5
            a(L)=1;
        else  a(L)=0;
        end
    end

    a1=modulation(a,index);
    a2=reshape(a1,Tx_n,frame_length/index/Tx_n);
    %sigma = 1/sqrt(2*SNR(i)*index);%standard deviation 
    sigma = 1/sqrt(2*SNR(i)*index);%standard deviation
    AWGN_noise = sigma*(randn(Rx_n,frame_length/index/Tx_n)+j*randn(Rx_n,frame_length/index/Tx_n));
    %dec0 = zeros(Tx_n,frame_length/index/Tx_n);    %for ML
    dec1 = zeros(frame_length/index , index);      %for ZF+OSIC
    dec2 = zeros(frame_length/index , index);       %for ZF
    dec3 = zeros(frame_length/index , index);       %for MMSE
    dec4 = zeros(frame_length/index , index);      % for MMSE+OSIC    
    dec5=zeros(1,frame_length/index);               %for QR
    dec6=zeros(1,frame_length);                     %for 
    out_mapping=mapping(index);
    for col_idx = 1:frame_length/index/Tx_n

        H=(randn(Rx_n,Tx_n)+j*randn(Rx_n,Tx_n))/sqrt(2);    %均值为0，方差为1的复高斯信道，包络服从瑞利分布
        H=sqrtm(Rt)*H*sqrtm(Rr);
        h=H;   %  for MMSE+OSIC
        r=H*a2(:,col_idx)+AWGN_noise(:,col_idx);    %垂直编码
        G=pinv(H);  %for ZF
        G1=inv(H'*H+sigma.^2*eye(Tx_n))*H';     %for MMSE
        [gk1 p0]=min(sum(abs(G1).^2,2));    %for MMSE+IC
        [gk k0]=min(sum(abs(G).^2,2));  %for ZF+OSIC
        r_ic = r;   %for ZF+OSIC 
        % QR
        [Q,R]=qr(H);
         y=Q'*r;    % ' 表示复共扼转置,对于酉矩阵，转置*本身=I
         G2=y(Tx_n)/R(Tx_n,Tx_n)-out_mapping;   %例如R =       1.8515            -0.6117 - 0.8030i
                                                %               0            -2.0710          
                                                %               0                  0          
                                                %               0                  0    
                                                %因为Rx_n>=Tx_n,所以R矩阵形似如上，所以
                                                %R(Rx_n,Tx_n)用R(Tx_n,Tx_n)代
                                                %替
                                                
    
         [Qk q]=min(abs(G2).^2);
        dec5(Tx_n+(col_idx-1)*Tx_n)=out_mapping(q);%是对第一列的最后一行的数据解码
         for m=Tx_n:-1:2    %是对第一列的前面几行数据解码
            G2=(y(m-1)-R(m-1,Tx_n:-1:m)*dec5((Tx_n:-1:m)+(col_idx-1)*Tx_n).')/R(m-1,m-1)-out_mapping;%.'表示非共轭阵列转置
            [Qk q]=min(abs(G2).^2);
            dec5(m-1+(col_idx-1)*Tx_n)=out_mapping(q);
         end
%**************************************************
         %          ML
       %count=0:1:2^Tx_n-1;
       %count1=zeros(Tx_n,2^Tx_n);
       %a4=zeros(Rx_n,2^Tx_n);
       %for  n1=1:2^Tx_n
            
       %      d=[-1,1];
       %       c1=de2bi(count(n1),Tx_n);
       %      c2=d(c1+1);
        %      count1(:,n1)=c2.';
        %     a4(:,n1)=r-H*count1(:,n1);
          
       %  end; % n1

        %     [gk k2]=min(sum(abs(a4).^2,1));    %1是求每一列的和,因为只有一列数据，就是求F范数
        %      dec0(:,col_idx)=count1(:,k2);
 %************************************************    
  
      
        %           % ZF
        det_zf = G*r;   %size(det_zf)= Tx_n  1
        dec2((col_idx-1)*Tx_n+[1:Tx_n],:) = reshape(demodulation(det_zf.',index),index,Tx_n).';
        %           MMSE
         det_MMSE=G1*r;
         dec3((col_idx-1)*Tx_n+[1:Tx_n],:) = reshape(demodulation(det_MMSE.',index),index,Tx_n).';
 
         
  
        %ZF+OSIC
        for m=1:Tx_n
           k1(m)=k0;
           y=G(k1(m),:)*r_ic;
           dec1((col_idx-1)*Tx_n+k1(m),:) = demodulation(y,index);      %解调后就可以排除了干扰，得到原始数据
           b=modulation(dec1((col_idx-1)*Tx_n+k1(m),:),index);          %然后再将原始数据映射到对应的调制后的数据，再在下面减去
          r_ic = r_ic - b*H(:, k1(m));
            H(:, k1(m))=zeros(Rx_n,1);		%接受分集体现在此?????????
            G=pinv(H);
            temp = sum(abs(G).^2,2);
            temp(k1(1:m)) = 1e10;   %e的10次方，inf
           [gk k0]=min(temp);
 
        end; %m
        
        %************************************************
         %MMSE+OSIC
          for n=1:Tx_n
            p1(n)=p0;
            y1=G1(p1(n),:)*r;
            dec4((col_idx-1)*Tx_n+p1(n),:) = demodulation(y1,index);
            b1=modulation(dec4((col_idx-1)*Tx_n+p1(n),:),index);
            r = r - b1*h(:, p1(n));
            h(:, p1(n))=zeros(Rx_n,1);
            G1=inv(h'*h+sigma.^2*eye(Tx_n))*h';
            temp1 = sum(abs(G1).^2,2);
            temp1(p1(1:n)) = 1e10;
            [gk1 p0]=min(temp1);
           end; %m
    end % col_idx
  %************************************************  
    %dec0=reshape(dec0,1,frame_length/index);
    %dec0=demodulation(dec0,index); 
    %NumErr0 =sum(dec0~=a)
    dec6=demodulation(dec5,index); 
    NumErr1 =sum(abs(reshape(dec1.',1,frame_length)~=a));
    NumErr2= sum(abs(reshape(dec2.',1,frame_length)~=a));
    NumErr3= sum(abs(reshape(dec3.',1,frame_length)~=a));
    NumErr4= sum(abs(reshape(dec4.',1,frame_length)~=a));
    NumErr5= sum(dec6~=a);
   % Ber0(i) = NumErr0/frame_length;
    Ber1(i) = NumErr1/frame_length;
    Ber2(i) = NumErr2/frame_length;
    Ber3(i) = NumErr3/frame_length;
    Ber4(i) = NumErr4/frame_length;
    Ber5(i) = NumErr5/frame_length;
    save ('data_4x4V_QPSK.mat', 'SNRindB','Ber1','Ber2','Ber3','Ber4','Ber5');
end % Eb/N0
%semilogy (SNRindB,Ber0,'<y-')
hold on
semilogy (SNRindB,Ber1,'*m-')
semilogy (SNRindB,Ber2,'+r-')
semilogy (SNRindB,Ber3,'og-')
semilogy (SNRindB,Ber4,'xk-')
semilogy (SNRindB,Ber5,'>c-')

xlabel('Eb/No in dB');
ylabel('error probability');
title('Tx = 4,Rx = 4,QPSK modulation');
legend('ZF+IS','ZF','MMSE','MMSE+IS','QR');
grid on