目录
零、前言
这是我研一通信系统仿真的专题答辩内容，就是当做笔记记录的，如果有内容上的错误请及时私信我，我会做出修改的，本文代码是可以用的，自己要多调试调试。
最近一段时间发现每天都有大量同学通过这篇博客加我问相关问题，我创建了一个交流群，专门讨论导频估计及通讯系统仿真等相关问题，想入群交流的关注同名公众号"通信汪的美好生活"或者扫描本文末尾二维码，然后私信我就好了。
一、为什么要信道估计
定义：用导频获取信道状态信息的这个过程就叫做信道估计。
从上图可知，与单纯的OFDM模型相比在里面的OFDM信道模型明显多出了插入导频和信道估计这两个过程，所以我们现在要解决的问题是为什么要进行信道估计。以及导频是干啥的。
信号在传输过程中会经历 1. 时间轴上的随机衰落（由于移动多普勒效应引起）见图2 2. 频率轴上的衰落（由于多径效应引起）。见图3
所以信道估计的目的，就是在于尽可能利用较少的，在时间方向和频率方向分布的pilot信号，对时间变化和频率变化进行抽样，然后利用这种对衰落已完备描述的抽样信息，经过滤波降噪和内插，恢复全频率（数据载波）的频率响应。信道估计基本的两种方法，一种是盲估计，一种是基于训练的估计，导频一般就应用在后者。
二、导频的概念
（1）为什么要有导频
在正式发送数据之前，先给系统训练训练，让它知道周围的通讯环境是怎样的，然后在正式发送数据时，发射端就可以根据这些信息(称为信道状态信息)做出相应的调整，从而达到较好的通信性能。
（2）导频在信道估计中作用
信道估计目前理解应该就是估计信道系数h，即如果发送信号已知，比如说x，理论上的接收信号应该是y=h*x+n，假设n=0，即不考虑噪声，那么信道系数就可以等于h=y/x。如果这个时候我们可以测量接收信号y^，那么就可以获得信道系数h^=y^/x。同时在一般使用中，我们强调导频都是一个序列，即一个已知序列，然后我们是采用相关的方法来进行上述的操作的。主要好处在于，首先相当于多个时刻的值做叠加，更稳定一些，同时由于噪声与设置的序列一般是无关的，那么其相关值即趋近与0，所以也能够避免噪声的干扰，增加准确性。
（3）关于导频序列的补充
在802.11协议中，基本都是固定子载波作为导频，并且用其作为信道估计。同时除了导频，802.11协议中也会采用LTF（长训练字段）做一个粗估计，方法也是同样采用固定的序列。 在LTE中，导频貌似和参考信号RS是一样的，尤其是小区专用参考信号CRS。 同时有关导频的内容在5G里面也比较被关注，其主要在于导频污染这个话题上。由于解调信号我们需要知道信道系数，同时信道系数就是利用导频所获得，如果导频互相干扰，那么就会严重通信质量。导频污染主要是由于不同的导频序列不正交导致（学术点可以认为序列间的汉明距离比较小），从而不同的序列做相关互相干扰，从而不能够很好估计信道。
三、最小二乘法估计
（1）LS信道估计算法分析
LS算法就是在不考虑噪声的条件下，估计信道的冲击响应向量hLS=[h0，h1，…，hN-1]T，使其代价函数最小。LS估计器的代价函数定义如下：
（2）LS信道估计的特点
LS信道估计的特点是简单，但是从其代价函数可以看出，在找最优解时没有考虑接收信号中的噪声及子载波间的干扰，因为这种算法估计出的信道对接收信号进行抽头复系数均衡时，输出信号的均方误差较大，准确度受到限制。
四、最小均方误差估计
（1）MMSE信道估计算法分析
MMSE信道估计算法对于ICI(子载波间干扰)和高斯白噪声有很好的抑制作用，它是在LS估计的基础上进行的：
（2）MMSE信道估计特点
MMSE算法的运算量要比LS算法大的多，随子载波数N呈指数增长，并需随导频信号X的变化实时进行矩阵的逆运算((XXH)一1)，导致系统效率很低。
（3）关于公式的推导
关于最小均方误差的公式推导，下面评论区链接提供了两种方法，感兴趣的可以自己去看。
五、简单运用
（一）MIMO_OFDM系统导频估计的MSE曲线图
请画出在信噪比从0-20dB下，MIMO_OFDM系统中随机导频序列设计的导频估计的MSE（Mean Square Error）曲线图。
要求：
1．MSE：h的估计值与h的期望值间差值的二范数，在循环多次下的均值（随着信噪比变化，求出误差，画出来）
2. OFDM块长（傅里叶变换长度）与各用户信道的多径数之和应存在什么关系？（ppt上回答，与最小二乘有关，看有没有唯一解）答：“各用户的多径数之和<=OFDM块长”才是正确答案，这涉及到的是未知数和方程个数之间的关系。
代码：
function Y=OFDM(x,h)
x_ifft=ifft(x);
x_cir=[x_ifft(length(x_ifft)-length(h):length(x_ifft));x_ifft;];
y_cir=conv(x_cir,h);
y_nocir=y_cir(length(h)+2:length(h)+2+length(x)-1);
Y=fft(y_nocir);clc
clear all;
close all;
NBlock=50;
WN=exp(-j*2*pi/NBlock); %DFT 离散傅里叶变换
for n=1:NBlock
for k=1:NBlock
DFTMatrix(k,n)=WN^((n-1)*(k-1)); %傅里叶变换矩阵构造
end
end
for i=1:7
X(:,i)=randn(NBlock,1);
end
%%用户数为7，并根据信道参数构造A
h1=[0.7;0.5;0.3;];
h2=[0.7;0.5;0.3;];
h3=[0.4;0.6;0.1;];
h4=[0.3;0.4;];
h5=[0.3;0.5;0.7;];
h6=[0.1;0.6;0.1;];
h7=[0.7;0.4;0.2;];
A=[diag(X(:,1))*DFTMatrix(:,1:length(h1)) diag(X(:,2))*DFTMatrix(:,1:length(h2)) diag(X(:,3))*DFTMatrix(:,1:length(h3)) diag(X(:,4))*DFTMatrix(:,1:length(h4)) diag(X(:,5))*DFTMatrix(:,1:length(h5)) diag(X(:,6))*DFTMatrix(:,1:length(h6)) diag(X(:,7))*DFTMatrix(:,1:length(h7))];
XMse(:,1)=zeros(NBlock,1)+1;
for K=1:NBlock
XMse(K,2)=exp(-j*2*pi*(length(h1))*(K-1)/NBlock);
end
for K=1:NBlock
XMse(K,3)=exp(-j*2*pi*(length(h1)+length(h2))*(K-1)/NBlock);
end
for K=1:NBlock
XMse(K,4)=exp(-j*2*pi*(length(h1)+length(h2)+length(h3))*(K-1)/NBlock);
end
for K=1:NBlock
XMse(K,5)=exp(-j*2*pi*(length(h1)+length(h2)+length(h3)+length(h4))*(K-1)/NBlock);
end
for K=1:NBlock
XMse(K,6)=exp(-j*2*pi*(length(h1)+length(h2)+length(h3)+length(h4)+length(h5))*(K-1)/NBlock);
end
for K=1:NBlock
XMse(K,7)=exp(-j*2*pi*(length(h1)+length(h2)+length(h3)+length(h4)+length(h5)+length(6))*(K-1)/NBlock);
end
AMse=[diag(XMse(:,1))*DFTMatrix(:,1:length(h1)) diag(XMse(:,2))*DFTMatrix(:,1:length(h2)) diag(XMse(:,3))*DFTMatrix(:,1:length(h3)) diag(XMse(:,4))*DFTMatrix(:,1:length(h4)) diag(XMse(:,5))*DFTMatrix(:,1:length(h5)) diag(XMse(:,6))*DFTMatrix(:,1:length(h6)) diag(XMse(:,7))*DFTMatrix(:,1:length(h7)) ];
%%随机导频
h=[h1;h2;h3;h4;h5;h6;h7];
for i=1:800
for snr=0:20
Y=A*h;
Y=awgn(Y,snr,'measured');
hEst=pinv(A)*Y; %信道估计
% LS_MSE(snr+1)= norm(hEst-mean(hEst*(1:length(hEst))))./21;
LS_MSE(snr+1,i)= norm((hEst)-h',2)/100;
end
end
LS_MSE800=mean( LS_MSE,2);
%%正交导频
for i=1:800
for snr=0:20
YMse=AMse*h;
YMse=awgn(YMse,snr,'measured');
hEstMse=pinv(AMse)*YMse;%信道估计h
%MMSE_MSE(snr+1)= norm(hEstMse-mean(hEstMse*(1:length(hEstMse))))./21;
MMSE_MSE(snr+1,i)= norm((hEstMse)-h',2)/100;
end
%MMSE_MSE()
end
MMSE_MSE800=mean( MMSE_MSE,2);
figure(1);
semilogy(0:20,LS_MSE800,'b-o',0:20,MMSE_MSE800,'r-*');
legend('随机导频','正交导频');
title('Mse与snr的关系图');
xlabel('snr');
ylabel('Mse');
（二）画出信噪比0：20下的误码率曲线
根据估计出的信道，将信号进行解调，画出信噪比0：20下的误码率曲线。（仅限于SISO_ofdm系统下的随机导频和正交导频）（注：在单输入单输出下）
要求：不能用if或case语句判决。
代码：
clc
clear all;
close all;
num=8000;
x_int=randi(8,num,1)-1; % 随机产生800个8进制整数
plot(x_int,'*');
for m=1:8
psk(m,1)=cos(11*pi/8-(m-1)*pi/4)+i*sin(11*pi/8-(m-1)*pi/4);
end
x=psk(x_int+1);% 8psk调制 x是x_int调制后的信号
BlockLength=8;
NBlock=8;
WN=exp(-j*2*pi/NBlock); %DFT 离散傅里叶变换
for n=1:NBlock
for k=1:NBlock
DFTMatrix(k,n)=WN^((n-1)*(k-1)); %傅里叶变换矩阵构造
end
end
h=[0.8;0.5;0.3;0.15;0.1;0.05];
BER_R=[];
%%随机导频信道估计
for i=1:100
for snr=0:30
X(:,1)=randn(NBlock,1);
A=[diag(X(:,1))*DFTMatrix(:,1:length(h))];
Y=A*h;
Y=awgn(Y,snr,'measured');
hEst=pinv(A)*Y; %信道估计h
y_awgn=awgn(x,snr,'measured');
y_multipath=conv(x,h);
for n=1:num/BlockLength
y_ofdm((n-1)*BlockLength+1:n*BlockLength,1)=OFDM(x((n-1)*BlockLength+1:n*BlockLength),h);
% 经过ofdm调制
end
y_ofdm=awgn(y_ofdm,snr,'measured');
HEst=fft([hEst;zeros(BlockLength-length(hEst),1)]);
for n=1:num/BlockLength
y_demod((n-1)*BlockLength+1:n*BlockLength)=y_ofdm((n-1)*BlockLength+1:n*BlockLength)./HEst;
end
%%系统判决
s=[0:BlockLength-1];
wholey_demod=[];
for mn=1:num
wholeD=[];
for m=1:BlockLength
D=norm((y_demod(mn)-psk(m)),2);
wholeD=[wholeD,D];
end
dom=min(wholeD);
w=find(wholeD==min(wholeD));
y_demod(mn)=s(w);
wholey_demod=[wholey_demod,y_demod(mn)];
end
wholey_demod1=wholey_demod';
symble_error=sum(not(wholey_demod1==x_int));
SER1=symble_error/num;
BER_R(i,snr+1)=SER1;%[BER_R,SER1];
end
end
BER_R800=mean(BER_R,1);
%%正交导频信道估计
BER_Q=[];
for i=1:100
for snr=0:30
XMse(:,1)=zeros(NBlock,1)+1;
AMse=[diag(XMse(:,1))*DFTMatrix(:,1:length(h))];
YMse=AMse*h;
YMse=awgn(YMse,snr,'measured');
hEstMse=pinv(AMse)*YMse;%信道估计h
y_awgn=awgn(x,snr,'measured');
y_multipath=conv(x,h);
for n=1:num/BlockLength
y_ofdm((n-1)*BlockLength+1:n*BlockLength,1)=OFDM(x((n-1)*BlockLength+1:n*BlockLength),h);
% 经过ofdm调制
end
y_ofdm=awgn(y_ofdm,snr,'measured');
HEstMse=fft([hEstMse;zeros(BlockLength-length(hEstMse),1)]);
for n=1:num/BlockLength
y_demod((n-1)*BlockLength+1:n*BlockLength)=y_ofdm((n-1)*BlockLength+1:n*BlockLength)./HEstMse;
end
%%系统判决
s=[0:BlockLength-1];
wholey_demod=[];
for mn=1:num
wholeD=[];
for m=1:BlockLength
D=norm((y_demod(mn)-psk(m)),2);
wholeD=[wholeD,D];
end
dom=min(wholeD);
w=find(wholeD==min(wholeD));
y_demod(mn)=s(w);
wholey_demod=[wholey_demod,y_demod(mn)];
end
wholey_demod1=wholey_demod';
symble_error=sum(not(wholey_demod1==x_int));
SER1=symble_error/num;
BER_Q(i,snr+1)=SER1;
end
end
BER_Q800=mean(BER_Q,1);
snr=0:30;
figure(1);
semilogy(snr,BER_R800,'b-o',snr,BER_Q800,'r-*');
legend('随机导频','正交导频');
title('BER与snr的关系图');
xlabel('snr');
ylabel('BER');六、结果分析
（1）为啥要进行蒙特卡洛？
通过大量的重复实验把某一件事情发生的频率当作概率。 如果重复次数少得到的结论是非常不可靠的。
（2）五里面题目第一问：
（3）五里面题目第二问：
七、结尾
利用MATLAB程序仿真实现了LS以及MMSE两种分别基于随机导频和正交导频的信道估计方法， 给出了两者的MSE及BER随信噪比变化的曲线。通过上面图片分析效果很明显，也不用多解释了，有什么疑问可以加我好友。我是一名热爱生活的通信汪，虽然生活很苦但是知识很甜，一起加油吧各位通信汪！！！
扫一扫
1003
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
