（内容源自详解MATLAB／SIMULINK 通信系统建模与仿真 刘学勇编著第九章内容，有兴趣的读者请阅读原书）
clear all
%%%%%%%参数设计部分%%%%%%%
Nsp=52;%系统子载波数（不包括直流载波）
Nfft=64;%FFT长度
Ncp=16;%循环前缀长度
Ns=Nfft+Ncp;%一个完整OFDM长度
noc=53;%包含直流载波的总的子载波数
Nd=6;%每帧包含的OFDM符号数（不包括训练符号）
M1=4;%QPSK调制
M2=16;%16QAM调制
sr=250000;%OFDM符号速率
EbN0=0:2:30;%归一化信噪比
Nfrm=10000;%每种信噪比下的仿真帧数
ts=1/sr/Ns;%OFDM抽样时间间隔
t=0:ts:(Ns*(Nd+1)*Nfrm-1)*ts;%抽样时刻
fd=100;%最大多普勒平移
h=rayleigh(fd,t);%生成单径瑞利信道
%训练符号频域数据，采用802.11a中的长训练符号数据
Preamble=[1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1];
Preamble1=zeros(1,Nfft);
Preamble1(2:27)=Preamble(27:end);%训练符号重拍后的数据
Preamble1(39:end)=Preamble(1:26);
preamble1=ifft(Preamble1);%训练符号时域数据
preamble1=[preamble1(Nfft-Ncp+1:end) preamble1];%加入循环前缀
首先来看参数设计部分，
载波的话一共有53个，其中一共有1个直流载波+52个系统子载波，其中52个系统子载波又分为48个数据子载波和4个导频子载波（见表9-1的2 3 4行）
FFT长度=64（见表9-1FFT周期TFFT）
循环前缀长度（见表9-1）
完整的OFDM长度（见表9-1，OFDM符号间隔）
每帧包含的OFDM符号数（见图9-14 PPDU帧结构，可以看到一个帧=2+1+？个OFDM符号，这里数据段data应该有3个OFDM符号，所以2+1+3=6，题目中也明确指出了6个符号构成一帧）
OFDM符号速率（见表9-1，1/OFDM符号间隔）
OFDM抽样时间间隔（1/sr/Ns=1/250000/80=50*10^-9，与表9-1抽样间隔一致）
抽样时刻：首先1个ts就是表9-1中的一个chip，一帧里面有6个OFDM数据符号和1个长训练符号（Nd+1），一个OFDM符号长度为80chip（见表9-1），一共的仿真帧数为10000(Nfrm)帧
所以最后一共含有的chip数为：Ns*(Nd+1)*Nfrm=N
所以t=0:ts:(N-1)*ts;
单径瑞利衰落信道：该信道的建模与时间t和最大多普勒频移fd有关。
训练符号频域数据：采用802.11a中的长符号训练数据（式9-16）
这里所使用的训练符号长度为64，其中前26个符号使用802.11a中的长训练符号数据的后26位，后26个符号使用数据的前26位(重排)，中间的全部用0进行填充。
最后加入的循环前缀是训练符号的后16位，在加入循环前缀之后，训练符号的长度为80，和一个OFDM符号的长度一致。
------------------------------------------------------------------------------------------------------------------------------
%%%%%%%仿真循环%%%%%%%
for ii=1:length(EbN0)
%******发射机部分******
msg1=randsrc(Nsp,Nd*Nfrm,[0:M1-1]);%QPSK信息数据
msg2=randsrc(Nsp,Nd*Nfrm,[0:M2-1]);%16-QAM信息数据
data1=pskmod(msg1,M1,pi/4);%QPSK调制
data2=qammod(msg2,M2)/sqrt(10);%16QAM调制并归一化
data3=zeros(Nfft,Nd*Nfrm);%根据FFT要求，对数据重排
data4=zeros(Nfft,Nd*Nfrm);
data3(2:27,:)=data1(27:end,:);
data3(39:end,:)=data1(1:26,:);
data4(2:27,:)=data2(27:end,:);
data4(39:end,:)=data2(1:26,:);
clear data1 data2;%清除不必要的临时变量
data3=ifft(data3);%IFFT变换
data4=ifft(data4);
data3=[data3(Nfft-Ncp+1:end,:);data3];%加入循环前缀
data4=[data4(Nfft-Ncp+1:end,:);data4];
spow1=norm(data3,'fro').^2/(Nsp*Nd*Nfrm);%计算数据符号能量
spow2=norm(data4,'fro').^2/(Nsp*Nd*Nfrm);
data5=zeros(Ns,(Nd+1)*Nfrm);%加入训练符号
data6=data5;
for indx=1:Nfrm
data5(:,(indx-1)*(Nd+1)+1)=preamble1.';
data5(:,(indx-1)*(Nd+1)+2 :indx*(Nd+1))=data3(:,(indx-1)*Nd+1:indx*Nd);
data6(:,(indx-1)*(Nd+1)+1)=preamble1.';
data6(:,(indx-1)*(Nd+1)+2 :indx*(Nd+1))=data4(:,(indx-1)*Nd+1:indx*Nd);
end
clear data3 data4
data5=reshape(data5,1,Ns*(Nd+1)*Nfrm);%并串变换
data6=reshape(data6,1,Ns*(Nd+1)*Nfrm);
sigma1=sqrt(1/2*spow1/log2(M1)*10.^(-EbN0(ii)/10));%根据EbN0计算噪声标准差
sigma2=sqrt(1/2*spow2/log2(M2)*10.^(-EbN0(ii)/10));
for indx=1:Nfrm
dd1=data5((indx-1)*Ns*(Nd+1)+1:indx*Ns*(Nd+1));%当前帧的发射数据(9)
dd2=data6((indx-1)*Ns*(Nd+1)+1:indx*Ns*(Nd+1));
hh=h((indx-1)*Ns*(Nd+1)+1:indx*Ns*(Nd+1));%当前帧对应的信道参数（10）
%信号通过单径瑞利衰落信道，并加入高斯白噪声（11）
r1=hh.*dd1+sigma1*(randn(1,length(dd1))+j*randn(1,length(dd1)));
r2=hh.*dd2+sigma2*(randn(1,length(dd2))+j*randn(1,length(dd2)));
r1=reshape(r1,Ns,Nd+1);%串并变换（12）
r2=reshape(r2,Ns,Nd+1);
r1=r1(Ncp+1:end,:);%移除循环前缀（13）
r2=r2(Ncp+1:end,:);接下来是发射机部分：
首先利用随机函数产生信息数据，msg1，msg2
randsrc(m,n,alphabet)：产生m行n列的随机数，alphabet是矩阵中元素可能取值的集合，也是等概率分布，alphabet中的重复的元素要忽略不计。
这里产生的数据有Nsp行（每一行对应OFDM中的一个子载波，共52行，直流载波是不传输数据的），
具体到每一个单独的行，共有Nd*Nfrm个数据，也就是每帧包含的OFDM符号数乘上仿真帧数，
即每个子载波上传输的OFDM符号数。
接下来进行QPSK和16QAM调制（data1,data2）
QPSK初始相位为pi/4，16QAM的归一化处理见下图
这里引用了阿__星的文章内容，文章链接如下
接下来是根据FFT要求，对数据进行重排，
重排的原因可以见图9-22上面的文字，具体来说就是
将53个子载波（前）映射到64个子载波（后）之上，映射的规则为
之前的训练符号频域数据的重排也是同理，因为训练符号也需要进行FFT运算。
之后加上循环前缀
计算数据符号能量
这里使用的函数为norm，使用方法如下：
n= norm(X,"fro") 返回矩阵或数组 X 的 Frobenius 范数。
关于此类型范数更多的信息，可以访问链接
https://zhuanlan.zhihu.com/p/680400515
之后的Nsp*Nd*Nfrm其实就是data3/data4矩阵含有的元素总数，也就是对范数元素总数的均值。
加入训练符号（6）
总结就是在每一帧前面都插入一个训练符号
并串变换（7）
根据EbN0计算噪声标准差（8）
（没看懂，待补）
当前帧的发送数据（9）
这里每次循环都取其中的一帧，因为之前已经进行了串并变换，所以现在的data5/6是一个向量，选择其中的一帧就需要从向量中选择Ns*(Nd+1)的长度
当前帧对应的信道参数（10）
这里我们使用向量t建立了一个长度为Ns*(Nd+1)*Nfrm的瑞利衰落信道向量h，可以发现h的长度与data5/6的长度是一致的，所以两者向量元素一一对应的关系，我们对其中的一帧经过信道时，也采用用信道中对应的部分
通过信道，并添加高斯白噪声（11）
这里通过信道采用的是点乘的方式，也就是说是信道中的元素与信号中的元素依次对应相乘，
例：信号【1 2 3 4 5】，信道【5 4 3 2 1】，经过信道的信号为【 5 8 9 8 5】，
换句话说，这里的信道起到的是改变信号幅值的作用。
这里在添加噪声的时候添加的是复高斯白噪声，因为QPSK和QAM都是复数信号。
串并变换（12）
这里将一帧信号由向量形式转成矩阵形式，其中矩阵的长为一帧中含有的符号个数（7），矩阵的宽为80（一个完整的OFDM长度）。
移除循环前缀（13）
之前在每一帧信号之前都添加了循环前缀，现在将循环前缀去掉（每一帧的前16行，只留下信号）
%%%%%%理想信道估计%%%%%%
hh=reshape(hh,Ns,Nd+1);%信道参数数据重排（1）
hh=hh(Ncp+1:end,:);
x1=r1(:,2:end)./hh(:,2:end);%信道补偿（2）
x2=r2(:,2:end)./hh(:,2:end);
x1=fft(x1);%fft运算
x2=fft(x2);
x1=[x1(39:end,:);x1(2:27,:)];%数据重排
x2=[x2(39:end,:);x2(2:27,:)];
x1=pskdemod(x1,M1,pi/4);%数据解调
x2=qamdemod(x2*sqrt(10),M2);
%统计一帧中的错误比特数
[neb1(indx),temp]=biterr(x1,msg1(:,(indx-1)*Nd+1:indx*Nd),log2(M1));
[neb2(indx),temp]=biterr(x2,msg2(:,(indx-1)*Nd+1:indx*Nd),log2(M2));
%%%%%%根据训练符号进行的信道估计%%%%%%%
R1=fft(r1);%fft运算
R2=fft(r2);
R1=[R1(39:end,:);R1(2:27,:)];%数据重排
R2=[R2(39:end,:);R2(2:27,:)];
HH1=(Preamble.')./R1(:,1);%信道估计
HH2=(Preamble.')./R2(:,1);
HH1=HH1*ones(1,Nd);%根据信道估计结果进行信道补偿
HH2=HH2*ones(1,Nd);
x3=R1(:,2:end).*HH1;
x4=R2(:,2:end).*HH2;
x3=pskdemod(x3,M1,pi/4);%数据解调
x4=qamdemod(x4.*sqrt(10),M2);
%统计一帧中的错误比特数
[neb3(indx),temp]=biterr(x3,msg1(:,(indx-1)*Nd+1:indx*Nd),log2(M1));
[neb4(indx),temp]=biterr(x4,msg2(:,(indx-1)*Nd+1:indx*Nd),log2(M2));
end
ber1(ii)=sum(neb1)/(Nsp*log2(M1)*Nd*Nfrm);%理想信道估计的误比特率
ber2(ii)=sum(neb2)/(Nsp*log2(M2)*Nd*Nfrm);
ber3(ii)=sum(neb3)/(Nsp*log2(M1)*Nd*Nfrm);%根据训练符号信道估计的误比特率
ber4(ii)=sum(neb4)/(Nsp*log2(M2)*Nd*Nfrm);
end接下来是理想信道估计部分
（1）信道参数数据重排，这里将信道也变换成一帧数据的矩阵形式
（2）信道补偿，信道补偿是指在通信系统中，通过一定的技术手段对信号在传输过程中受到的信道影响进行补偿，以减少信号失真和提高通信质量。
例：Y=XH，则X=Y/H
注意，这里没有对第一列进行信道补偿，因为第一列是训练符号，之后的才是信息，是需要补偿的。
（3）fft运算，数据重排（数据重排过程中将第一列的训练符号舍弃了）
（4）数据解调（这里qam在调制的时候采用了归一化，这里消去归一化的影响）
（5）统计一帧中的错误比特数(待补)
接下来是根据训练符号进行的信道估计
这里可以看到，训练符号进行的信道估计和理想信道估计的区别在于
理想信道估计是直接把信道矩阵拿过来用（将信道定为已知量）
训练信道进行的信道估计是利用训练信号间接的得到信道矩阵的情况，信道是未知量
用已知的训练信号点除经过信道后的训练信道，可以得到信道情况(HH1/2)
HH1=HH1*ones(1,Nd)构建信道矩阵
利用信道矩阵进行信道补偿
结合以上几部分的完整源码见下：
clear all
%%%%%%%参数设计部分%%%%%%%
Nsp=52;%系统子载波数（不包括直流载波）
Nfft=64;%FFT长度
Ncp=16;%循环前缀长度
Ns=Nfft+Ncp;%一个完整OFDM长度
noc=53;%包含直流载波的总的子载波数
Nd=6;%每帧包含的OFDM符号数（不包括训练符号）
M1=4;%QPSK调制
M2=16;%16QAM调制
sr=250000;%OFDM符号速率
EbN0=0:2:30;%归一化信噪比
Nfrm=10000;%每种信噪比下的仿真帧数
ts=1/sr/Ns;%OFDM抽样时间间隔
t=0:ts:(Ns*(Nd+1)*Nfrm-1)*ts;%抽样时刻
fd=100;%最大多普勒平移
h=rayleigh(fd,t);%生成单径瑞利信道
%训练符号频域数据，采用802.11a中的长训练符号数据
Preamble=[1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1];
Preamble1=zeros(1,Nfft);
Preamble1(2:27)=Preamble(27:end);%训练符号重拍后的数据
Preamble1(39:end)=Preamble(1:26);
preamble1=ifft(Preamble1);%训练符号时域数据
preamble1=[preamble1(Nfft-Ncp+1:end) preamble1];%加入循环前缀
%%%%%%%仿真循环%%%%%%%
for ii=1:length(EbN0)
%******发射机部分******
msg1=randsrc(Nsp,Nd*Nfrm,[0:M1-1]);%QPSK信息数据
msg2=randsrc(Nsp,Nd*Nfrm,[0:M2-1]);%16-QAM信息数据
data1=pskmod(msg1,M1,pi/4);%QPSK调制
data2=qammod(msg2,M2)/sqrt(10);%16QAM调制并归一化
data3=zeros(Nfft,Nd*Nfrm);%根据FFT要求，对数据重排
data4=zeros(Nfft,Nd*Nfrm);
data3(2:27,:)=data1(27:end,:);
data3(39:end,:)=data1(1:26,:);
data4(2:27,:)=data2(27:end,:);
data4(39:end,:)=data2(1:26,:);
clear data1 data2;%清除不必要的临时变量
data3=ifft(data3);%IFFT变换
data4=ifft(data4);
data3=[data3(Nfft-Ncp+1:end,:);data3];%加入循环前缀
data4=[data4(Nfft-Ncp+1:end,:);data4];
spow1=norm(data3,'fro').^2/(Nsp*Nd*Nfrm);%计算数据符号能量
spow2=norm(data4,'fro').^2/(Nsp*Nd*Nfrm);
data5=zeros(Ns,(Nd+1)*Nfrm);%加入训练符号
data6=data5;
for indx=1:Nfrm
data5(:,(indx-1)*(Nd+1)+1)=preamble1.';
data5(:,(indx-1)*(Nd+1)+2 :indx*(Nd+1))=data3(:,(indx-1)*Nd+1:indx*Nd);
data6(:,(indx-1)*(Nd+1)+1)=preamble1.';
data6(:,(indx-1)*(Nd+1)+2 :indx*(Nd+1))=data4(:,(indx-1)*Nd+1:indx*Nd);
end
clear data3 data4
data5=reshape(data5,1,Ns*(Nd+1)*Nfrm);%并串变换
data6=reshape(data6,1,Ns*(Nd+1)*Nfrm);
sigma1=sqrt(1/2*spow1/log2(M1)*10.^(-EbN0(ii)/10));%根据EbN0计算噪声标准差
sigma2=sqrt(1/2*spow2/log2(M2)*10.^(-EbN0(ii)/10));
for indx=1:Nfrm
dd1=data5((indx-1)*Ns*(Nd+1)+1:indx*Ns*(Nd+1));%当前帧的发射数据
dd2=data6((indx-1)*Ns*(Nd+1)+1:indx*Ns*(Nd+1));
hh=h((indx-1)*Ns*(Nd+1)+1:indx*Ns*(Nd+1));%当前帧对应的信道参数
%信号通过单径瑞利衰落信道，并加入高斯白噪声
r1=hh.*dd1+sigma1*(randn(1,length(dd1))+j*randn(1,length(dd1)));
r2=hh.*dd2+sigma2*(randn(1,length(dd2))+j*randn(1,length(dd2)));
r1=reshape(r1,Ns,Nd+1);%串并变换
r2=reshape(r2,Ns,Nd+1);
r1=r1(Ncp+1:end,:);%移除循环前缀
r2=r2(Ncp+1:end,:);
%%%%%%理想信道估计%%%%%%
hh=reshape(hh,Ns,Nd+1);%信道参数数据重排
hh=hh(Ncp+1:end,:);
x1=r1(:,2:end)./hh(:,2:end);%信道补偿
x2=r2(:,2:end)./hh(:,2:end);
x1=fft(x1);%fft运算
x2=fft(x2);
x1=[x1(39:end,:);x1(2:27,:)];%数据重排
x2=[x2(39:end,:);x2(2:27,:)];
x1=pskdemod(x1,M1,pi/4);%数据解调
x2=qamdemod(x2*sqrt(10),M2);
%统计一帧中的错误比特数
[neb1(indx),temp]=biterr(x1,msg1(:,(indx-1)*Nd+1:indx*Nd),log2(M1));
[neb2(indx),temp]=biterr(x2,msg2(:,(indx-1)*Nd+1:indx*Nd),log2(M2));
%%%%%%根据训练符号进行的信道估计%%%%%%%
R1=fft(r1);%fft运算
R2=fft(r2);
R1=[R1(39:end,:);R1(2:27,:)];%数据重排
R2=[R2(39:end,:);R2(2:27,:)];
HH1=(Preamble.')./R1(:,1);%信道估计
HH2=(Preamble.')./R2(:,1);
HH1=HH1*ones(1,Nd);%根据信道估计结果进行信道补偿
HH2=HH2*ones(1,Nd);
x3=R1(:,2:end).*HH1;
x4=R2(:,2:end).*HH2;
x3=pskdemod(x3,M1,pi/4);%数据解调
x4=qamdemod(x4.*sqrt(10),M2);
%统计一帧中的错误比特数
[neb3(indx),temp]=biterr(x3,msg1(:,(indx-1)*Nd+1:indx*Nd),log2(M1));
[neb4(indx),temp]=biterr(x4,msg2(:,(indx-1)*Nd+1:indx*Nd),log2(M2));
end
ber1(ii)=sum(neb1)/(Nsp*log2(M1)*Nd*Nfrm);%理想信道估计的误比特率
ber2(ii)=sum(neb2)/(Nsp*log2(M2)*Nd*Nfrm);
ber3(ii)=sum(neb3)/(Nsp*log2(M1)*Nd*Nfrm);%根据训练符号信道估计的误比特率
ber4(ii)=sum(neb4)/(Nsp*log2(M2)*Nd*Nfrm);
end
semilogy(EbN0,ber1,'-ro',EbN0,ber3,'-rv',EbN0,ber2,'-r*',EbN0,ber4,'-rd')
grid on
title('OFDM系统误比特率性能')
legend('QPSK理想信道估计','QPSK训练符号信道估计','16-QAM理想信道估计','16-QAM训练符号信道估计')
xlabel('信噪比(EbN0)')
ylabel('误比特率')最后绘制图像
结论为：
P.S.:这里是简化版的仿真，没有考虑扰码和卷积编码等内容
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
