stu_vblast.m

VBLST译码算法仿真源程序 代码包括了常用的调制和解调方式。

close all;
clear all;
clc;
Tx_n = 2;               %4x4 天线配置
index = 4;              %调制方式:表示为1/2/4/6 比特每符号
frame_length = 48;      %帧长
numloop=1000;           %循环次数
SNR_dB_MIN=0;           %信噪比区间
SNR_dB_MAX=14;
SNR_dB_stmp=2;          %信噪比步进
bertmp=zeros(1,numloop);%预定义临时比特差错矩阵
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;               %最终的比特差错矢量
for num=1:length(Rx_num);
    Rx_n=Rx_num(num);
if Rx_n < Tx_n          %如果接收天线数多于发送天线，报错
    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                %外循环，以信噪比步进为间隔
    for looptmp=1:numloop                                   %内循环，同一信噪比下，循环numloop次
        source = rand(1,Tx_n*frame_length*index)>.5;        %产生信源序列，1行4x96x2比特，这里没有考虑信道编码了
        column_num = length(source) / Tx_n;                 %VBLAST分层，每层的比特数
        encoder_output = zeros(Tx_n, column_num);           %内编码输出是4x(96*2)的矩阵
        for col = 1:column_num                              %VBLAST比特分配——内编码
            for row = 1:Tx_n
                encoder_output(row, col) = source((col-1)*Tx_n + row);
            end
        end
        modu_output=zeros(Tx_n,frame_length);               %调制输出，符号矩阵为4x96
        for t=1:Tx_n                                        %调用调制函数
            modu_output(t,:) = modulation(encoder_output(t,:), index);
        end
        H=zeros(Rx_n,Tx_n);                                 %信道矩阵预定义
        for rx = 1:Rx_n
            for tx = 1:Tx_n
                randn('state',sum(100*clock));              %每次重置到不同状态
                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));                        %信噪比对应的噪声方差
        r = decoder_input;                                  %译码输入就是接收信号
        y = zeros(Tx_n,frame_length);

        G1=pinv(H);                                         %求伪逆,ZF
        H1=H;                                               %保留一次信道矩阵，用于ZF
        r1=r;                                               %保留一次接收信号，用于ZF
        G=(H'*H+delta2*eye(Tx_n))\H';                       %MMSE对应的G矩阵
        [gk k0]=min(sum(abs(G).^2,2));                      %该运算确定检测的顺序
        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);     %解调矩阵预定义
        demodu_output1 =demodu_output ;
        for t = 1:Tx_n                                      %调用解调函数
            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分层空时解码
            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);                      %计算本轮循环的比特差错
        [Num1,Ber1]=symerr(sink1,source);
        bertmp(looptmp)=Ber;                                %存进临时比特差错矩阵
        bertmp1(looptmp)=Ber1;
    end
    BER_mmse(num,(SNR_dB-SNR_dB_MIN)/SNR_dB_stmp+1)=mean(bertmp);    %对应每个信噪比，平均numloop
    BER_zf(num,(SNR_dB-SNR_dB_MIN)/SNR_dB_stmp+1)=mean(bertmp1);  %次循环的比特差错结果
end
end

figure;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_mmse(1,:),'r*:');      %2发8收
hold on;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_mmse(2,:),'g*:');      %4发8收
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_mmse(3,:),'*:');      %8发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发8收
hold on;
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_zf(2,:),'g*:');      %4发8收
semilogy(SNR_dB_MIN:SNR_dB_stmp:SNR_dB_MAX,BER_zf(3,:),'*:');      %8发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');