% 論文仿真原型程序,仿真四種（ZF,ZF-SIC,MMSE,MMSE-SIC） Vblast接收機(jī)的檢測性能，繪制誤比特率～信噪比曲線。
% Zhaidewei 2007.4.18 (2007-5-20修正版)
% Email:zeal0502@yahoo.com.cn
% qq:402731924

% 發(fā)端初始化===============================================================
% 發(fā)射天線數(shù)tx,接收天線數(shù)rx,發(fā)射矩陣長度L(幀長)
tx=4;rx=4;L=10000;
Modulation='BPSK';
EbNo=[0:5:20];
B=30000;Ts=1/24300;
% 建立EbN0與SNR之間的換算關(guān)系
SNR=EbNo-10*log10(Ts*B);
% 信源A
A=randint(tx*L,1);
% 經(jīng)過BPSK調(diào)制的V-Blast發(fā)射矩陣X
X=zeros(tx,L);
for k=1:tx
X(k,:)=((-1).^(A(k:tx:end)+1))';
end

% 信道傳輸=================================================================

% 快衰落Rayleigh信道H
H=sqrt(1/2)*(randn(rx,tx,L)+i*randn(rx,tx,L));
% 均值為0方差為1的高斯白噪聲n
n=sqrt(1/2)*(randn(rx,L)+i*randn(rx,L));
% 未疊加噪聲的接收信號R
R=zeros(rx,L);
for k=1:L
    R(:,k)=H(:,:,k)*X(:,k);
end

% 檢測

%ZF========================================================================
disp('berz');
berz=[];
% 在不同的信噪比下計(jì)算ZF接收機(jī)誤比特率berz
  for m=SNR
      m
      % 每個(gè)子信道的平均信噪比為snr的接受信號R_noised
      snr=10^(m/10);
      R_noised=R+sqrt(1/snr)*n;
      x=[];
      a=zeros(tx*L,1);
      % 逐時(shí)隙對接收符號矢量進(jìn)行檢測，合并得到一幀發(fā)射矩陣X的估計(jì)x
      for t=1:L
          r=R_noised(:,t);
          % 迫零矩陣G
          G=pinv(H(:,:,t));
          y=G*r;
          xtemp=(y>=0)-(y<0)+0;
          x=[x,xtemp];
      end
      
      % 從x求A的估計(jì)a
      for k=1:tx
          a(k:tx:end)=(x(k:tx:end)+1)/2;
      end
      % 比較A和a計(jì)算錯(cuò)值率temp_ber
      [errbit,temp_ber]=biterr(A,a);
      berz=[berz,temp_ber];
  end
  figure
  semilogy(EbNo,berz,'o- b')
 
% ZF-SIC(ordered)==========================================================
disp('ZF-SIC');
berzs=[];
  for m=SNR
      m
      snr=10^(m/10);
      R_noised=R+sqrt(1/snr)*n;
      x=[];
      a=zeros(tx*L,1);
      for t=1:L
          r=R_noised(:,t);
          HH=H(:,:,t);
          G=pinv(HH);
          S=[1:tx];% S表示一個(gè)時(shí)隙內(nèi)還未檢測的符號的序號的集合
          xtemp=zeros(tx,1);
          % 逐發(fā)射天線進(jìn)行檢測
          for k=1:tx
              % G的在序號集合S的非零元素集合內(nèi)范數(shù)最小的行是wki，它是G的第ki行
              [wki,ki]=minnorm(G,S);
              % 已經(jīng)檢測過的序號清零
              S(ki)=0;
              % 判決統(tǒng)計(jì)量y
              y=wki*r;
              % BPSK判決
              xtemp(ki)=1*(y>=0)-1*(y<0)+0;
              % SIC串行干擾抵消
              r=r-xtemp(ki)*H(:,ki,t);
              % 將已經(jīng)檢測的信號對應(yīng)的信道矩陣的列清零
              HH(:,ki)=zeros(rx,1);
              G=pinv(HH);
          end
          x=[x,xtemp];
      end
      for k=1:tx
          a(k:tx:end)=(x(k:tx:end)+1)/2;
      end
      [errbit,temp_ber]=biterr(A,a);
      berzs=[berzs,temp_ber];
  end
  hold on
  semilogy(EbNo,berzs,'o- r')

 % MMSE====================================================================
 disp('MMSE')
  berm=[];
  for m=SNR
      m
      snr=10^(m/10);
      R_noised=R+sqrt(1/snr)*n;
      x=[];
      a=zeros(tx*L,1);
      for t=1:L
          r=R_noised(:,t);
          HH=H(:,:,t);
          xtemp=zeros(tx,1);
          w=inv(HH'*HH+(1/snr)*eye(tx))*HH';
          y=w*r;
          xtemp=(y>=0)-(y<0)+0;
          x=[x,xtemp];
      end
      for k=1:tx
          a(k:tx:end)=(x(k:tx:end)+1)/2;
      end
      [errbit,temp_ber]=biterr(A,a);
      berm=[berm,temp_ber];
  end
  semilogy(EbNo,berm,'s- b')
 % MMSE-SIC(ordered)=======================================================
 disp('MMSE-SIC')
 berms=[];
  for m=SNR
      m
      snr=10^(m/10);
      R_noised=R+sqrt(1/snr)*n;
      x=[];
      a=zeros(tx*L,1);
      for t=1:L
          r=R_noised(:,t);
          HH=H(:,:,t);
          G=pinv(HH);
          S=[1:tx];
          xtemp=zeros(tx,1);
          for k=1:tx
              [wki,ki]=minnorm(G);
              w=inv(HH'*HH+(1/snr)*eye(tx-k+1))*HH';
              y=w(ki,:)*r;
              xtemp(S(ki))=(y>=0)-(y<0)+0;
              r=r-xtemp(S(ki))*H(:,S(ki),t);
              HH(:,ki)=[];
              S(ki)=[];
              G=pinv(HH);
          end
          x=[x,xtemp];
      end
      for k=1:tx
          a(k:tx:end)=(x(k:tx:end)+1)/2;
      end
      [errbit,temp_ber]=biterr(A,a);
      berms=[berms,temp_ber];
  end
  semilogy(EbNo,berms,'s- m'),grid on
 legend('ZF','ZF-SIC','MMSE','MMSE-SIC');
 xlabel('Eb/No(dB)');
 ylabel('誤比特率');
 title('blast檢測比較')