核心提示：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
微信扫一扫：
