基于导频的信道估计实现-阿里云开发者社区

零、前言

这是我研一通信系统仿真的专题答辩内容，就是当做笔记记录的，如果有内容上的错误请及时私信我，我会做出修改的，本文代码是可以用的，自己要多调试调试。

一、为什么要信道估计

定义：用导频获取信道状态信息的这个过程就叫做信道估计。

微信图片_20221017100318.png

从上图可知，与单纯的OFDM模型相比在里面的OFDM信道模型明显多出了插入导频和信道估计这两个过程，所以我们现在要解决的问题是为什么要进行信道估计。以及导频是干啥的。

信号在传输过程中会经历 1. 时间轴上的随机衰落（由于移动多普勒效应引起）见图2 2. 频率轴上的衰落（由于多径效应引起）。见图3

微信图片_20221017100323.png

图2移动多普勒效应

微信图片_20221017100325.png

图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信道估计的特点是简单，但是从其代价函数可以看出，在找最优解时没有考虑接收信号中的噪声及子载波间的干扰，因为这种算法估计出的信道对接收信号进行抽头复系数均衡时，输出信号的均方误差较大，准确度受到限制。

微信图片_20221017100332.png

四、最小均方误差估计

（1）MMSE信道估计算法分析

MMSE信道估计算法对于ICI(子载波间干扰)和高斯白噪声有很好的抑制作用，它是在LS估计的基础上进行的：

微信图片_20221017100335.png

微信图片_20221017100339.png

（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块长”才是正确答案，这涉及到的是未知数和方程个数之间的关系。

代码：

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）五里面题目第一问：

微信图片_20221017101059.png

信噪比0：20，用户数7，多径数20，蒙特卡洛100次

（3）五里面题目第二问：

微信图片_20221017101102.png

信噪比0：20，用户数1，多径数6，蒙特卡洛100次

七、结尾

利用MATLAB程序仿真实现了LS以及MMSE两种分别基于随机导频和正交导频的信道估计方法，给出了两者的MSE及BER随信噪比变化的曲线。通过上面图片分析效果很明显，也不用多解释了，有什么疑问可以加我好友。我是一名热爱生活的通信汪，虽然生活很苦但是知识很甜，一起加油吧各位通信汪！！！

基于导频的信道估计实现