sphereandml.m

word文档对球形译码的理论原理和仿真结果进行了分析 1.sphereandML是主程序

close all;
clear all;

%/////////////////////////初始值///////////////////////////////////%
snr=0:3:15;                                              %snr的dB形式  
snr_1=zeros(1,length(snr));                              %用来存储snr正常形式
sumber=0;                                                %初始化sumber
sumber1=0;                                               %初始化sumber1
sumber2=0;                                               %初始化sumber2
SymTime=10^3;                                            %仿真次数
mod_scheme='4qam';                                       %调制方式
mod_num=4;                                               %调制方式参数
avrgber=zeros(1,length(snr));                            %初始化误比特率，2范数球形译码
avrgber1=zeros(1,length(snr));                           %初始化误比特率，无穷范数球形译码
avrgber2=zeros(1,length(snr));                           %初始化误比特率，ML 
Nr=4;                                                    %接收天线数目
Nt=4;                                                    %发送天线数目
Number=10;                                               %bit source 
t1=zeros(1,length(snr));                                 %初始化记录时间t1
t2=zeros(1,length(snr));                                 %初始化记录时间t2
t3=zeros(1,length(snr));

for n=1:SymTime
    H1(:,:,n)=(randn(Nr,Nt)+j*randn(Nr,Nt))/sqrt(2);
end                                                      %产生信道矩阵H  

for ISNR=1:length(snr)
    snr_1(ISNR)=10^(snr(ISNR)/10);                       %将dB转化为正常形式                                         
    sumber=0;                                            %记录ber的和
    sumber1=0;
    sumber2=0;
    for round=1:SymTime                                      
        msg=randint(log2(mod_num),Nt*Number/2);
        h=modem.qammod('M',mod_num,'InputType','bit','symbolorder','Gray');%qam调制
        s=modulate(h,msg);
        constelPoints=[1+j 1-j -1+j -1-j];                      %$ for ML detect
  %     constelPoints=[-3-3*j -3-j -3+3*j -3+j -1-3*j -1-j -1+3*j -1+j 3-3*j 3-j 3+3*j 3+j 1-3*j 1-j 1+3*j 1+j]; 
        s1=[real(s);imag(s)];                                    %s1为信号的实矩阵
        s2=reshape(s,Nt,Number/2);                               %s2为信号的Nt行的复矩阵
        if mod_num==4
           codebook=[-1 1];                                       %表示在变为实矩阵后，4qam只有两种情况1和-1 
        end
        if mod_num==16                                          %表示在变为实矩阵后，16qam有4种情况
           codebook=[-3 -1 1 3]; 
        end   
        ss=reshape(s1,2*Nt,Number/2);                            %ss为信号的2*Nt行的实矩阵  
%       H=[real(H1(:,:,round)) -imag(H1(:,:,round));imag(H1(:,:,round)) real(H1(:,:,round))];                %将信道复矩阵变为信道实矩阵 


%/////////////////////////////////////////////////////////////////%
       r1=H1(:,:,round)*s2;                                 
       for ii=1:Nt                                              %添加噪声
           rr(ii,:)=awgn(r1(ii,:),snr_1(ISNR),'measured');
       end                                                      %rr为接收复信号 
     
       C1=3*(2*Nt)*2*(2/(10^(snr(ISNR)/10)));                                   %初始半径 二范数4qam
       C2=16*2*(2/(10^(snr(ISNR)/10)));                                      %初始半径无穷范数
%        C1=3*(2*Nt)*(10/(10^(snr(ISNR)/10)));                               %初始半径    二范数16qam  
%        C2=8*2*(10/(10^(snr(ISNR)/10)));                                      %初始半径  无穷范数
       rev1=[real(rr);imag(rr)];                                      %rr转变为实信号
       H=[real(H1(:,:,round)) -imag(H1(:,:,round));imag(H1(:,:,round)) real(H1(:,:,round))];              %信道实矩阵
     
     
       for tmp=1:Number/2                                            %二范数球形译码
           rev=rev1(:,tmp);   
           tic
           y_2norm_temp = spheredecode(rev, C1, H, codebook);
           for jj=1:Nt
               y_2norm(Nt*tmp+jj-Nt)=y_2norm_temp(jj)+j*y_2norm_temp(jj+Nt);                 %存储为复数形式
           end
           t1(ISNR)=toc+t1(ISNR); 
       end   
     
       H=[real(H1(:,:,round)) -imag(H1(:,:,round));imag(H1(:,:,round)) real(H1(:,:,round))];              %信道实矩阵
       for tmp=1:Number/2                                            %二范数球形译码
           rev=rev1(:,tmp);   
           tic
           y_infnorm_temp = spheredecodeinf(rev, C2, H, codebook);    
           for jj=1:Nt
               y_infnorm(Nt*tmp+jj-Nt)= y_infnorm_temp(jj)+j* y_infnorm_temp(jj+Nt);                 %存储为复数形式
           end
           t2(ISNR)=toc+t2(ISNR); 
       end   
     
     
     
       for tmp=1:Number/2
           min=100000;
           tic
           for flag1=1:mod_num
               for flag2=1:mod_num
                   for flag3=1:mod_num
                       for flag4=1:mod_num
                           stmp=[constelPoints(flag1);constelPoints(flag2);constelPoints(flag3);constelPoints(flag4)];
                           rtmp=rr(:,tmp);
%                   stmp=[constelPoints(flag1);constelPoints(flag2)];
%                   rtmp=[rr(1,tmp);rr(2,tmp)];
                           Maximum=norm(rtmp-H1(:,:,round)*stmp); 
                           if Maximum<min
                              min=Maximum;
                              slast=stmp;
                           end    
                       end
                   end
               end
           end  
           for jj=1:Nt        
               y_ML(Nt*tmp+jj-Nt)=slast(jj);
           end     
           t3(ISNR)=toc+t3(ISNR); 
        end  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
       h=modem.qamdemod('M',mod_num,'OutputType','bit','symbolorder','Gray');   %第二范数
       DM1_msg=demodulate(h,y_2norm);
       DM_msg=reshape(DM1_msg,log2(mod_num),Nt*Number/2);
       [number,ratio] = biterr(msg,DM_msg);  

       h=modem.qamdemod('M',mod_num,'OutputType','bit','symbolorder','Gray');
       DM3_msg=demodulate(h,y_infnorm);
       DM2_msg=reshape(DM3_msg,log2(mod_num),Nt*Number/2);
       [number1,ratio1] = biterr(msg,DM2_msg);                           %无穷范数  

       h=modem.qamdemod('M',mod_num,'OutputType','bit','symbolorder','Gray');
       DM5_msg=demodulate(h,y_ML);
       DM4_msg=reshape(DM5_msg,log2(mod_num),Nt*Number/2);
       [number2,ratio2] = biterr(msg,DM4_msg);                           %ML 
   
       sumber=sumber+ratio;  
       sumber1=sumber1+ratio1;
       sumber2=sumber2+ratio2;   
    end %round=1：SymTime
    avrgber(ISNR)=sumber/SymTime;
    avrgber1(ISNR)=sumber1/SymTime;
    avrgber2(ISNR)=sumber2/SymTime;
end

%******************plot*********************************************
figure(1)
semilogy(snr,avrgber,'-^b',snr,avrgber1,'-*g',snr,avrgber2,'-or')                %蓝色是2norm球形译码，绿色是无穷范数球形译码，红色是ML检测
title('average BER')
xlabel('SNR(dB)')
ylabel('BER')
grid on
figure(2)
plot(snr,t1,'-^b',snr,t2,'-*g',snr,t3,'-or')                %蓝色是2norm球形译码，绿色是无穷范数球形译码，红色是ML检测
title('simulation time')
xlabel('SNR(dB)')
ylabel('TIME')
grid on