clc;
clear;
close all;
len_s=100000; % 信号长度
snr=0:2:20;%信噪比
len_snr=length(snr);
s2_1=zeros(1,4);
ber_zf=zeros(1,len_snr);
ber_zf_sic=zeros(1,len_snr);
ber_zf_chsic=zeros(1,len_snr);
for ii=1:len_snr
    error_zf2=0;
  for i=1:len_s
    sn=randi([0,1],4,1); %产生随机信号
    mod_sn=pskmod(sn,2);
    sigma =sqrt(0.5.*10.^(-snr(ii)./10));
    %I,Q两路独立产生噪声
    Q_noise = randn(4, 1);
    I_noise = randn(4, 1);
    noise = complex(I_noise,Q_noise);
    noise = sigma.*noise;%高斯白噪声
    H=1/sqrt(2)*(randn(4,4) + 1i*randn(4,4));
    y1=H*mod_sn+noise;
    G=inv(H);  %求逆  
    y1_1=G*y1;
    s=pskdemod(y1_1,2);%BPSK解调
    error_zf1=find(s~=sn);%错误判决
    error_zf2=length(error_zf1)+error_zf2;
  end
    ber_zf(ii)=error_zf2./(len_s*4);
end

for ii=1:len_snr
    error_zf22=0;
    for i=1:len_s
       sn2=randi([0,1],4,1); %产生随机信号
       mod_sn2=pskmod(sn2,2);
       sigma =sqrt(0.5.*10.^(-snr(ii)./10));
       %I,Q两路独立产生噪声
       Q_noise = randn(4, 1);
       I_noise = randn(4, 1);
       noise = complex(I_noise,Q_noise);
       noise = sigma.*noise;%高斯白噪声
       H=1/sqrt(2)*(randn(4,4) + 1i*randn(4,4));
       y2_1_1=H*mod_sn2+noise;  %4*1
       for p=1:4
          G2=pinv(H); %伪逆矩阵
          w2_1=G2(p,:);  %1*4
          y2_1_2=w2_1*y2_1_1;
          s2_1(p)=pskdemod(y2_1_2,2);%BPSK解调
          s2_1_1=pskmod(s2_1(p),2);   %先检测第一个信号
          aH=s2_1_1.*H(:,p);
          y2_1_1= y2_1_1-aH;
          H(:,p)=zeros(4,1);
       end
      error_zf21=find(s2_1.'~=sn2);%错误判决
      error_zf22=length(error_zf21)+error_zf22;
    end
    ber_zf_sic(ii)=error_zf22./(len_s*4);
end
%选信号强度最大的
for ii=1:len_snr
    error_zf22=0;
    for i=1:len_s
       sn2=randi([0,1],4,1); %产生随机信号
       mod_sn2=pskmod(sn2,2);  %bpsk调制
       sigma =sqrt(0.5.*10.^(-snr(ii)./10));
       %I,Q两路独立产生噪声
       Q_noise = randn(4, 1);
       I_noise = randn(4, 1);
       noise = complex(I_noise,Q_noise);
       noise = sigma.*noise;%高斯白噪声
       H=1/sqrt(2)*(randn(4,4) + 1i*randn(4,4));
       y2_1_1=H*mod_sn2+noise;
       G3_1=zeros(4,4);
       for p=1:4
          G3=pinv(H);%伪逆矩阵
          sum1=sum(abs(G3+G3_1).^2,2);%对矩阵行求和
          [x,p1]=min(sum1); %
          G3_1(p1,:)=[1000 1000 1000 1000];
          w2_1=G3(p1,:);
          y2_1_2=w2_1*y2_1_1;
          s2_1(p1)=pskdemod(y2_1_2,2);%BPSK解调
          s2_1_1=pskmod(s2_1(p1),2);
          aH=s2_1_1.*H(:,p1);
          y2_1_1= y2_1_1-aH;
          H(:,p1)=zeros(4,1);
       end
      error_zf21=find(s2_1.'~=sn2);%错误判决
      error_zf22=length(error_zf21)+error_zf22;
    end
    ber_zf_chsic(ii)=error_zf22./(len_s*4);
end
 semilogy(snr,ber_zf,'-k*',snr,ber_zf_sic,'-ko',snr,ber_zf_chsic,'-kx')
 grid on
 legend('迫零检测算法','串行干扰消除算法','排序串行干扰消除算法')
 xlabel('信噪比snr')
 ylabel('误比特率')