close all;
clear all;
clc;
Tx_n = 2;               %4x4 天線配置
index = 4;              %調(diào)制方式:表示為1/2/4/6 比特每符號(hào)
frame_length = 48;      %幀長(zhǎng)
numloop=1000;           %循環(huán)次數(shù)
SNR_dB_MIN=0;           %信噪比區(qū)間
SNR_dB_MAX=14;
SNR_dB_stmp=2;          %信噪比步進(jìn)
bertmp=zeros(1,numloop);%預(yù)定義臨時(shí)比特差錯(cuò)矩陣
bertmp1=bertmp;
Rx_num=[2 4 8];
BER_mmse=zeros(length(Rx_num),length(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX));
BER_zf=BER_mmse;               %最終的比特差錯(cuò)矢量
for num=1:length(Rx_num);
    Rx_n=Rx_num(num);
if Rx_n < Tx_n          %如果接收天線數(shù)多于發(fā)送天線，報(bào)錯(cuò)
    disp_str = ('error!!!the Rx_n must no less than Tx_n');
    disp(disp_str);
    return;
end
for SNR_dB=SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX                %外循環(huán)，以信噪比步進(jìn)為間隔
    for looptmp=1:numloop                                   %內(nèi)循環(huán)，同一信噪比下，循環(huán)numloop次
        source = rand(1,Tx_n*frame_length*index)>.5;        %產(chǎn)生信源序列，1行4x96x2比特，這里沒有考慮信道編碼了
        column_num = length(source) / Tx_n;                 %VBLAST分層，每層的比特?cái)?shù)
        encoder_output = zeros(Tx_n, column_num);           %內(nèi)編碼輸出是4x(96*2)的矩陣
        for col = 1:column_num                              %VBLAST比特分配——內(nèi)編碼
            for row = 1:Tx_n
                encoder_output(row, col) = source((col-1)*Tx_n + row);
            end
        end
        modu_output=zeros(Tx_n,frame_length);               %調(diào)制輸出，符號(hào)矩陣為4x96
        for t=1:Tx_n                                        %調(diào)用調(diào)制函數(shù)
            modu_output(t,:) = modulation(encoder_output(t,:), index);
        end
        H=zeros(Rx_n,Tx_n);                                 %信道矩陣預(yù)定義
        for rx = 1:Rx_n
            for tx = 1:Tx_n
                randn('state',sum(100*clock));              %每次重置到不同狀態(tài)
                x = randn(1);
                y = randn(1);
                alpha = sqrt(x^2+y^2);                      %Rayleigh分布
                theta = 2*pi*rand(1);
                H(rx,tx) = alpha*exp(j*theta);
            end
        end
        channel_output = H*modu_output;                     %信道輸出，H*s，還未加噪聲
        decoder_input = awgn(channel_output,SNR_dB,'measured');%y=a*H+n，加噪聲，作為收端譯碼輸入
        delta2 = 1/(10^(SNR_dB/10));                        %信噪比對(duì)應(yīng)的噪聲方差
        r = decoder_input;                                  %譯碼輸入就是接收信號(hào)
        y = zeros(Tx_n,frame_length);

        G1=pinv(H);                                         %求偽逆,ZF
        H1=H;                                               %保留一次信道矩陣，用于ZF
        r1=r;                                               %保留一次接收信號(hào)，用于ZF
        G=(H'*H+delta2*eye(Tx_n))\H';                       %MMSE對(duì)應(yīng)的G矩陣
        [gk k0]=min(sum(abs(G).^2,2));                      %該運(yùn)算確定檢測(cè)的順序
        k1=zeros(1,Tx_n);
        w=zeros(Tx_n,Rx_n);
        a=zeros(Tx_n,frame_length);
        for m=1:Tx_n                                        %基于MMSE的干擾抵消算法
            k1(m)=k0;
            w(m,:)=G(k1(m),:);
            y=w(m,:)*r;
            a(k1(m),:)=Q(y,index);
            r = r -H(:, k1(m))*a(k1(m),:);
            H(:,k0)=zeros(Rx_n,1);
            G=pinv(H);
            for t=1:m
                G(k1(t),:)=inf;
            end
            [gk k0]=min(sum(abs(G).^2,2));
        end
        decoder_output=a;

        [gk k0]=min(sum(abs(G1).^2,2)); %?????????????
        k1=zeros(1,Tx_n);
        w=zeros(Tx_n,Rx_n);
        a=zeros(Tx_n,frame_length);
        for m=1:Tx_n                                        %基于ZF的干擾抵消算法
            k1(m)=k0;
            w(m,:)=G1(k1(m),:);
            y=w(m,:)*r1;
            a(k1(m),:)=Q(y,index);
            r1 = r1 -H1(:, k1(m))*a(k1(m),:);
            H1(:,k0)=zeros(Rx_n,1);
            G1=pinv(H1);
            for t=1:m
                G1(k1(t),:)=inf;
            end
            [gk k0]=min(sum(abs(G1).^2,2));
        end
        decoder_output1=a;

        demodu_output = zeros(Tx_n,frame_length*index);     %解調(diào)矩陣預(yù)定義
        demodu_output1 =demodu_output ;
        for t = 1:Tx_n                                      %調(diào)用解調(diào)函數(shù)
            demodu_output(t,:) = demodulation(decoder_output(t,:), index);
            demodu_output1(t,:) = demodulation(decoder_output1(t,:), index);
        end
        sink = zeros(1,Tx_n*frame_length*index);            %信宿序列
        sink1 = sink;
        for i = 1:Tx_n                                      %VBLAST分層空時(shí)解碼
            sink(i:Tx_n:end)=demodu_output(i,:);
            sink1(i:Tx_n:end)=demodu_output1(i,:);
        end
        source=double(source);
        [Num,Ber]=symerr(sink,source);                      %計(jì)算本輪循環(huán)的比特差錯(cuò)
        [Num1,Ber1]=symerr(sink1,source);
        bertmp(looptmp)=Ber;                                %存進(jìn)臨時(shí)比特差錯(cuò)矩陣
        bertmp1(looptmp)=Ber1;
    end
    BER_mmse(num,(SNR_dB-SNR_dB_MIN)/SNR_dB_stmp+1)=mean(bertmp);    %對(duì)應(yīng)每個(gè)信噪比，平均numloop
    BER_zf(num,(SNR_dB-SNR_dB_MIN)/SNR_dB_stmp+1)=mean(bertmp1);  %次循環(huán)的比特差錯(cuò)結(jié)果
end
end

figure;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_mmse(1,:),'r*:');      %2發(fā)8收
hold on;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_mmse(2,:),'g*:');      %4發(fā)8收
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_mmse(3,:),'*:');      %8發(fā)8收
grid on;
xlabel('SNR(db)');ylabel('BER');
title('V-BLAST of MMSE');
legend('Tn=2,Rn=2','Tn=2,Rn=4','Tn=2,Rn=8');
figure(2);
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_zf(1,:),'r*:');      %2發(fā)8收
hold on;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_zf(2,:),'g*:');      %4發(fā)8收
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_zf(3,:),'*:');      %8發(fā)8收
grid on;
xlabel('SNR(db)');ylabel('BER');
title('V-BLAST of ZF');
legend('Tn=2,Rn=2','Tn=2,Rn=4','Tn=2,Rn=8');