您现在的位置:首页 >> ★免费资源 >> 源码下载 >> 内容

hslogic算法仿真-基于QAM调制的MIMO系统仿真

时间:2019-7-17 23:00:17 点击:

  核心提示:hslogic算法仿真-基于QAM调制的MIMO系统仿真...
clear;t0=clock;
 
% 发射天线数tx,接收天线数rx,发射矩阵长度L(帧长),一帧(所有发射天线)所包含的总比特数
tx=4; rx=4; %tx,rx可修改
Modulation='16QAM';
bpc=4*tx; L=bpc*50000;
%xAll是QAM信号的全集;aAll是QAM解调信号组合的全集
GrayTable=[4 3 1 2 8 7 5 6 16 15 13 14 12 11 9 10];
xAll=[-3-3i,-3-1i,-3+1i,-3+3i,   -1-3i,-1-1i,-1+1i,-1+3i,   1-3i,1-1i,1+1i,1+3i,   3-3i,3-1i,3+1i,3+3i ].';
aAll=[0 0 1 0; 0 0 1 1; 0 0 0 1; 0 0 0 0;
         0 1 1 0; 0 1 1 1; 0 1 0 1; 0 1 0 0;
         1 1 1 0; 1 1 1 1; 1 1 0 1; 1 1 0 0;
         1 0 1 0; 1 0 1 1; 1 0 0 1; 1 0 0 0];
EbNo=5:5:25;%dB
B=30000;%Hz
Rb=24300;%bits per second
% 建立EbN0与SNR之间的换算关系
SNR=EbNo+10*log10(Rb/B);
% 信源A
%A=randint(L,1);
A=randi([0,1],[L,1]);%R2010b and later version

% 经过16QAM调制的V-Blast发射矩阵X
X=zeros(tx,L/bpc);
b=zeros(L/bpc,1);
gray=zeros(L/bpc,1);

for k=1:tx
    b(:)=8*A((4*(k-1)+1):bpc:end)+4*A((4*(k-1)+2):bpc:end)+2*A((4*(k-1)+3):bpc:end)+A((4*(k-1)+4):bpc:end)+1;
    gray(:)=GrayTable(b(:));
    X(k,:)=xAll(gray(:));
end

% 信道传输=================================================================
% 快衰落Rayleigh信道H
H=sqrt(1/2)*(randn(rx,tx,L/bpc)+1i*randn(rx,tx,L/bpc));
%均值为0,方差为2.1180(与16QAM信号功率相同)的高斯白噪声n
N=(1+sqrt(5)/2)*(randn(rx,L/bpc)+1i*randn(rx,L/bpc));
% 未叠加噪声的接收信号R
R=zeros(rx,L/bpc);
for k=1:L/bpc
    R(:,k)=H(:,:,k)*X(:,k);
end

%QRM-MLD==================================================================
%Qt,Rt是Ht进行QR分解后得到的两个矩阵
%qdi是象限检测的次数,对于16QAM,需要进行3次象限检测
%qd(qdi)是每次象限检测变换坐标的偏移量
%qdr是象限检测划分的星座空间子区域,对于16QAM,共有2^3*2^3=64个子区域
%SRC是按欧氏距离排序的候选码元集合(Symbol Replica Candidate)
disp('QRM-MLD with ASESS')
Sm=[16,16,16,16];%Sm(m)是第m阶段的预设幸存候选码元个数

BER=[];
for msnr=SNR
    snr=10^(msnr/10);
    Rns=R+sqrt(1/snr)*N;

x=zeros(L/4,1);
a=zeros(L,1);
for t=1:L/bpc
t;
    Ht=H(:,:,t);     
    Nt=N(:,t);
    r=Rns(:,t);

    %================排序 
    HtExt=zeros(4+4,4);
    HtExt(1:4,:)=Ht;
    HtExt(5:8,:)=(1/snr)*eye(4);
    Rt=zeros(4,4); Qt=HtExt; Pt=eye(4,4);

    for ii=1:4
        for qi=ii:4
         qt(qi)=norm(Qt(:,qi));
        end
        [~,ki]=min(qt);
    
        exchangeQt=Qt(:,ii);exchangeRt=Rt(:,ii);exchangePt=Pt(:,ii);
        Qt(:,ii)=Qt(:,ki);Rt(:,ii)=Rt(:,ki);Pt(:,ii)=Pt(:,ki);
        Qt(:,ki)=exchangeQt;Rt(:,ki)=exchangeRt;Pt(:,ki)=exchangePt;
    
        Rt(ii,ii)=(norm(Qt(:,ii)));
        Qt(:,ii)=Qt(:,ii)/Rt(ii,ii);
        if (ii+1<=4)
            for l=ii+1:4
                Rt(ii,l)=Qt(:,ii)'*Qt(:,l);
                Qt(:,l)=Qt(:,l)-Rt(ii,l)*Qt(:,ii);
            end
        end
        qt(ii)=65535;%qt(ii)已用过,求最小时不再包括qt(ii)
    end
%================
%     [Qt,Rt]=qr(Ht);
%   

rExt=zeros(tx+rx,1);
rExt(1:rx)=r;

Qt=Qt(1:tx,:);
   Z=Qt'*r;
   M=ones(4,round(max(Sm)));
   Mi=ones(4,round(max(Sm)));
   e=ones(4,round(max(Sm)),round(max(Sm)));
   c1=ones(4,round(max(Sm)));
   c2=ones(4,round(max(Sm)));
   c3=ones(4,round(max(Sm)));
   c4=ones(4,round(max(Sm)));
   E=ones(4,round(max(Sm)));

    
%第一阶段
   temp(4)=Z(4)/Rt(4,4);
   SRC=quadrantdect(temp(4));
   Ed=Z(4)-Rt(4,4)*SRC(1:Sm(1));%未经计算的欧氏距离
%    Ed=Z(4)/Rt(4,4)-SRC(1:Sm(1));%未经计算的欧氏距离
   e(1,1:Sm(1),1)=real(Ed(1:Sm(1))).^2+imag(Ed(1:Sm(1))).^2;
   %计算欧氏距离e。e的第1个下标表示第几阶段,第2个下标表示本阶段幸存候选码元序号,第3个下标表示前一阶段幸存候选码元的序号.
   E(1,1:Sm(1))=e(1,1:Sm(1),1);
   c1(4,:)=SRC(1:Sm(1));
   M(1,1:Sm(1))=E(1,1:Sm(1));
  
   %第2-4阶段
       m=2;
       M(m,1:Sm(m-1))=M(m-1,1:Sm(m-1));
       Mi(m,1:Sm(m-1))=1;
       
       for j=1:1:Sm(m)
           [~,ci]=min(M(m,:));
           c2(4,j)=c1(4,ci);
           temp(3)=(Z(3)-Rt(3,4)*c2(4,j))/Rt(3,3);
           SRC=quadrantdect(temp(3));
           c2(3,j)=SRC(Mi(m,ci));
           dstc=(Z(3)-Rt(3,4)*c2(4,j))  -  Rt(3,3)*c2(3,j);
%            dstc=(Z(3)-Rt(3,4)*c2(4,j))/ Rt(3,3)- c2(3,j);
          
           e(2,ci,Mi(m,ci))=E(m-1,ci)+real(dstc).^2+imag(dstc).^2;
           E(m,j)=e(2,ci,Mi(m,ci));
           M(m,j)=E(m,j);
           Mi(m,j)=Mi(m,j)+1;
           
       end

       m=3;
       M(m,1:Sm(m-1))=M(m-1,1:Sm(m-1));
       Mi(m,1:Sm(m-1))=1;
       
       for j=1:Sm(m)
           [~,ci]=min(M(m,:));
           c3(4,j)=c2(4,ci);
           c3(3,j)=c2(3,ci);
%            c3(1:m-1,j)=c2(:,ci)
           temp(2)=(Z(2)-Rt(2,4)*c2(4,j)-Rt(2,3)*c2(3,j))/Rt(2,2);
           SRC=quadrantdect(temp(2));
           c3(2,j)=SRC(Mi(m,ci));
           dstc=(Z(2)-Rt(2,4)*c3(4,j)-Rt(2,3)*c3(3,j))  -  Rt(2,2)*c3(2,j);
%            dstc=(Z(2)-Rt(2,4)*c3(4,j)-Rt(2,3)*c3(3,j))/Rt(2,2)  -  c3(2,j);
           e(3,ci,Mi(m,ci))=E(m-1,ci)+real(dstc).^2+imag(dstc).^2;
           E(m,j)=e(3,ci,Mi(m,ci));
%            M(m,ci)=E(m,j);
%            Mi(m,ci)=Mi(m,ci)+1;
           M(m,j)=E(m,j);
           Mi(m,j)=Mi(m,j)+1;
           
       end

       m=4;
       
       j=1;
       M(m,1:Sm(m-1))=M(m-1,1:Sm(m-1));
       Mi(m,1:Sm(m-1))=1;
        [MM,ci]=min(M(m,:));
            c4(4,j)=c3(4,ci);
            c4(3,j)=c3(3,ci);
            c4(2,j)=c3(2,ci);
           temp(1)=(Z(1)-Rt(1,4)*c4(4,j)-Rt(1,3)*c4(3,j)-Rt(1,2)*c4(2,j))/Rt(1,1);
           SRC=quadrantdect(temp(1));
           c4(1,j)=SRC(Mi(m,ci));
           
%        for j=1:Sm(m)
%             [MM,ci]=min(M(m,:));
%             c4(4,j)=c3(4,ci);
%             c4(3,j)=c3(3,ci);
%             c4(2,j)=c3(2,ci);
% %            c4(1:m-1,j)=c3(:,ci)
%            temp(1)=(Z(1)-Rt(1,4)*c4(4,j)-Rt(1,3)*c4(3,j)-Rt(1,2)*c4(2,j))/Rt(1,1);
%            SRC=quadrantdect(temp(1));
%            c4(1,j)=SRC(Mi(m,ci));
%            dstc=(Z(1)-Rt(1,4)*c4(4,j)-Rt(1,3)*c4(3,j)-Rt(1,2)*c4(2,j))  -  Rt(1,1)*c4(1,j);
%            e(4,ci,Mi(m,ci))=E(m-1,ci)+real(dstc).^2+imag(dstc).^2;
%            E(m,j)=e(4,ci,Mi(m,ci));
% %            M(m,ci)=E(m,j);
% %            Mi(m,ci)=Mi(m,ci)+1;
%             M(m,j)=E(m,j);
%             Mi(m,j)=Mi(m,j)+1;
%            
%        end
c4=Pt*c4;%SQRD恢复顺序

x( tx*t-3 : tx*t )=c4(:,1);
end

ni=0;
 for idxS=1:L/4
     for idxT=1:16
         if x(idxS,1)==xAll(idxT,1)
             ni=idxT;
         end
    end
     a(4*(idxS-1)+1:4*idxS)=aAll(ni,:);
 end
     [number,bertemp] = biterr(A,a);
     BER=[BER,bertemp];
end
hold on
semilogy(SNR,BER,'>- b')


%ZF_SIC===================================================================
disp('ZF-SIC')
BER=[];
for m=SNR
    snr=10^(m/10);
    Rns=R+sqrt(1/snr)*N;
x=zeros(tx,L/bpc);
a=zeros(L,1);
    for t=1:L/bpc
          r=Rns(:,t);
          Ht=H(:,:,t);
          G=inv(Ht'*Ht)*Ht';
          S=(1:tx);% S表示第t个时隙内还未检测的符号的序号的集合
          xtemp=zeros(tx,1);
          atemp=zeros(1,4*tx);
          % 逐发射天线进行检测
          for k=1:tx
              %GG用来存储G的各行范数,以进行比较
              GG=zeros(size(G,1),1);
              for kk=1:size(G,1)  
                  GG(kk)=norm(G(kk,:));
              end
              %G的范数最小的行是wki,它是G的第ks行,第ks行对应的未检测信号序号为ki, 即S(ks)=ki
               [w,ks]=min(GG);
               ki=S(ks);
               wki=G(ks,:);
              % 将已经检测过的序号删除
              S(ks)=[];
              % 判决统计量y
              y=wki*r;
              %16QAM判决
              yv=ones(16,1)*y;%将y复制到一个16维列向量中
              dv=yv-xAll;
              dvn=abs(dv);
              [dni,ni]=min(dvn);
              xtemp(ki)=xAll(ni);
              atemp(4*(ki-1)+1:4*ki)=aAll(ni,:);
              % SIC串行干扰抵消
              r=r-xtemp(ki)*H(:,ki,t);
              % 将已经检测的信号对应的信道矩阵的列清零
              %Ht(:,ki)=zeros(rx,1);
              Ht(:,ks)=[];
              G=inv(Ht'*Ht)*Ht';
          end
          x(:,t)=xtemp;
          a(4*tx*(t-1)+1:4*tx*(t-1)+bpc)=atemp;
    end
    m
    [number,bertemp] = biterr(A,a);
    BER=[BER,bertemp];
end
%作图
semilogy(SNR,BER,'*- r')


% MMSE_SIC===================================================================
disp('MMSE-SIC')
BER=[];
for m=SNR
    snr=10^(m/10);
    Rns=R+sqrt(1/snr)*N;
x=zeros(tx,L/bpc);
a=zeros(L,1);
    for t=1:L/bpc
          r=Rns(:,t);
          Ht=H(:,:,t);
          G=inv(Ht'*Ht+(1/snr)*eye(tx))*Ht';
          S=(1:tx);% S表示第t个时隙内还未检测的符号的序号的集合
          xtemp=zeros(tx,1);
          atemp=zeros(1,4*tx);
          % 逐发射天线进行检测
          for k=1:tx
              %GG用来存储G的各行范数,以进行比较
              GG=zeros(size(G,1),1);
              for kk=1:size(G,1)  
                  GG(kk)=norm(G(kk,:));
              end
              %G的范数最小的行是wki,它是G的第ks行,第ks行对应的未检测信号序号为ki, 即S(ks)=ki
               [w,ks]=min(GG);
               ki=S(ks);
               wki=G(ks,:);
              % 将已经检测过的序号删除
              S(ks)=[];
              % 判决统计量y
              y=wki*r;
              %16QAM判决
              yv=ones(16,1)*y;%将y复制到一个16维列向量中
              dv=yv-xAll;
              dvn=abs(dv);
              [dni,ni]=min(dvn);
              xtemp(ki)=xAll(ni);
              atemp(4*(ki-1)+1:4*ki)=aAll(ni,:);
              % SIC串行干扰抵消
              r=r-xtemp(ki)*H(:,ki,t);
              % 将已经检测的信号对应的信道矩阵的列清零
              %Ht(:,ki)=zeros(rx,1);
              Ht(:,ks)=[];
              G=inv(Ht'*Ht+(1/snr)*eye(tx-k))*Ht';
          end
          x(:,t)=xtemp;
          a(4*tx*(t-1)+1:4*tx*(t-1)+bpc)=atemp;
    end
    m
    [number,bertemp] = biterr(A,a);
    BER=[BER,bertemp];
end
hold on
semilogy(SNR,BER,'+- g')
grid on 
%MLD=======================================================================
disp('MLD')

%生成一个时隙发射信号的全集,idxR=Row index idxC=Column index idxT=Table(xAll)index
C=zeros(tx,16^tx);
for idxR=1:tx
    for idxT=1:16
        C(idxR,16^(idxR-1)*(idxT-1)+1:16^(idxR):end)=xAll(idxT);
    end
end
for idxR=1:tx
    for idxC=1:16^tx
       if C(idxR,idxC)~=0
           c=C(idxR,idxC);
       else
           C(idxR,idxC)=c;
       end  
    end
end
%--------------------------------------------------------------------------
x=zeros(tx,L/bpc);
a=zeros(L,1);
xVct=zeros(1,L/4);

BER=[];
dC=zeros(rx,16^tx);
for m=SNR
    snr=10^(m/10);
    Rns=R+sqrt(1/snr)*N;

for t=1:L/bpc
    r=Rns(:,t);
    Ht=H(:,:,t);
    rC=Ht*C;
    avrgnC=zeros(1,16^tx);
    
    for idxC=1:rx
        dC(idxC,:)=rC(idxC,:)-r(idxC);
    end
    normC=sqrt(real(dC).^2+imag(dC).^2);
    for idxR=1:rx
        avrgnC=avrgnC+1/rx*normC(idxR,:);
    end
    [minid,idxR]=min(avrgnC);
    xtemp=C(:,idxR);
    x(:,t)=xtemp;
    xVct(tx*(t-1)+1:tx*t)=xtemp.';
end
for idxS=1:L/4
    for idxT=1:16
        if xVct(idxS)==xAll(idxT)
            ni=idxT;
        end
    end
    a(4*(idxS-1)+1:4*idxS)=aAll(ni,:);
end
m
    [number,bertemp] = biterr(A,a);
    BER=[BER,bertemp];
end
hold on
semilogy(SNR,BER,'>- b')
grid on

% legend('ZF','MMSE-SIC','MLD');
xlabel({'SNR(dB)'});
ylabel({'BER'});
% title('V-BLAST检测比较,16QAM调制')
elapsedTime=etime(clock,t0)


联系:highspeedlogic

QQ :1224848052

微信:HuangL1121

邮箱:1224848052@qq.com

网站:http://www.mat7lab.com/

网站:http://www.hslogic.com/

微信扫一扫:

Tags:QAM调制 
作者:QAM调制 来源:QAM调制
  • 您是如何找到本站的?
  • 百度搜索
  • Google搜索
  • 查阅资料过程中
  • 论坛发现
  • 百度贴吧发现
  • 朋友介绍
本站最新成功开发工程项目案例
相关评论
发表我的评论
  • 大名:
  • 内容:
  • matlab代做|matlab专业代做|matlab淘宝代做(www.hslogic.com) © 2019 版权所有 All Rights Reserved.
  • Email:highspeed_logic@163.com 站长QQ: 1224848052

    专业代做/代写/承接、MATLAB、SIMULINK、FPGA项目、博士/硕士/本科毕业设计、课题设计、论文,毕业论文,Coursework、Eassy、Assignment