vblast_cor.m

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

% 本程序是发射和接收天线之间均有相关性的情况下VBLAST系统
close all;
clear
echo off
Tx_n = 2;
Rx_n = 2;
index =2;
frame_length=1000;
SNRindB=0:2:20;
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.9;
                  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.9;
                  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); 
    AWGN_noise = sigma*(randn(Rx_n,frame_length/index/Tx_n)+j*randn(Rx_n,frame_length/index/Tx_n));
      dec1 = zeros(frame_length/index , index);
      dec2 = zeros(frame_length/index , index);
     dec3 = zeros(frame_length/index , index);
     dec4 = zeros(frame_length/index , index);
     
    for col_idx = 1:frame_length/index/Tx_n
       
        H=(randn(Rx_n,Tx_n)+j*randn(Rx_n,Tx_n))/sqrt(2);
          H1=sqrtm(Rr)*H*sqrtm(Rt);     %sqrtm求的是矩阵的平方根
         h=H1;
         r=H1*a2(:,col_idx)+AWGN_noise(:,col_idx);
          G=pinv(H1);
        G1=inv(H1'*H1+sigma.^2*eye(Tx_n))*H1';
         [gk1 p0]=min(sum(abs(G1).^2,2));
          [gk k0]=min(sum(abs(G).^2,2));
         r_ic = r;
        
      %  ZF
         det_zf = G*r;
         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*H1(:, k1(m));
            H1(:, k1(m))=zeros(Rx_n,1);
            G=pinv(H1);
            temp = sum(abs(G).^2,2);
            temp(k1(1:m)) = 1e10;
            [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
      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))
      Ber1(i) = NumErr1/frame_length;
      Ber2(i) = NumErr2/frame_length;
     Ber3(i) = NumErr3/frame_length;
     Ber4(i) = NumErr4/frame_length;
      save data_2x2_0.9CO_V_QPSK.mat SNRindB Ber1 Ber2 Ber3 Ber4;
end; % Eb/N0

% semilogy (SNRindB,Ber1,'-*')
% hold on
% semilogy (SNRindB,Ber2,'-+')
% semilogy (SNRindB,Ber3,'-o')
% semilogy (SNRindB,Ber4,'-.')


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