第7章：OFDM 信道估计与均衡（4）

本文所有可运行代码下载地址是：123kevin456/OFDM-

一、导频结构与图案

前三讲介绍了OFDM经过AWGN信道和衰落信道的误码率情况，其中在第（2）节假设了接收端已经完美知道信道状态信息h，做了完美补偿（perfect compensation）。

那接收端是如何知道信道状态信息的呢？知道了又有什么用呢？这就要从导频估计说起了。

首先回到熟悉的表达式：y=hx+n

假设只有一条径的衰落信道，不妨设x是bpsk调制，即可取1，或者-1。

发送端发送一大堆符号，为方便起见，比如 $x=\left[1,1,1,1\right]$ ，经过信道后，接收端收到 $y=\left[0.2,0.5,0.6,0.3\right]$ 。

此时一个问题问你，h等于多少呢？面对这一串数据，你会做什么样的处理呢？

可能你马上会想到 $\frac{sum(y)}{1+1+1+1}=0.4$ 。那么恭喜你，这便是最大似然估计（ML）。

（想想为什么是最大似然估计？以及说是最大似然估计，是否有前提条件要加上？思考几秒钟？）

当然，后面会介绍LS（最小平方）估计和MMSE（最小均方误差）估计。

那么导频符号是按照什么样的分布、出现在什么位置呢？这便是导频结构的问题。

导频结构可以分为三种：块状导频、梳状导频和格状导频。

看完上面的三种导频结构后，有两个关键问题出现在你的面前：

（1）导频数量：即需要插入多少个导频。

若插入导频数量多了，会消耗系统的时频资源，减少了有效数据的传输；插入导频少了，可能难以有效反应信道的真实特性；

（2）导频位置：即导频插在哪些地方。

比如是放置在第1号、5号、9号子载波位置呢？还是放置在2号、6号、10号子载波位置呢？有区别吗？

这当然有区别的哦。（体现在插值函数里，你可以看到区别）

本次仿真是采用梳状导频信号来做的信道估计，代码中可以看到导频的位置与间隔。

我设置的有16个导频符号，间隔为4。

二、基于导频的信道估计算法和插值方法

基于导频的信道估计方法，是在发送数据流中插入已知数据（导频符号），接收机根据这些已知数据经过信道衰落后的变化，与原始的导频符号进行比较，就得到了导频信号所在时刻和子载波上信道衰落的估计值。

而估计性能的优劣，一方面取决于使用算法的优劣，比如LS、MMSE算法，也取决于发送端设置的导频图案是否能够较好的反应信道的特性。

注意到上面这句“导频信号所在时刻和子载波上信道衰落的估计值”，那不是导频信号所在时刻和子载波信道，怎么办？

注意观察上面的导频结构图，均是时间和频率二维的。

不管是采用梳状、块状还是格状，都只知道其中某几个点的信道状况（对应到图里，便是黑色实心点），即要通过这些黑色点的信道状态信息去想办法知道白色点的信道状况。

这便要介绍插值了。

常见的插值函数是interp1，可以在MATLAB中help interp1中学习一下。下面给出的MATLAB代码也是使用的interp1这个插值函数。

除了上面基于导频的估计方法，其实ODFM还有基于判决反馈的信道估计方法、基于DFT的信道估计、基于叠加信号的信道估计、盲信道估计等等方法，在《MIMO-OFDM无线通信技术及MATLAB实现》中均有介绍，哈哈哈哈哈，这本书再一次被我点名表扬。

下面来看最重要的LS和MMSE均衡算法：

书上的6.14-6.17公式，我到没有想得非常明白。因为在6.14中，需要用到H，但接收端是不知道H的。

在一个指数衰减的多径功率时延谱（PDP）中，频域相关为6.16。那什么时候算是多径功率时延谱呢，如果不是指数衰减的多径功率时延谱，怎么进行处理呢？类似地问题，书中并没有讲清楚。

因此，在这部分代码中，我是直接调用书中自带的MMSE_CE.m这个代码的。

三、完整可运行MATLAB代码及其注意点

先展示实验结果：

图2 不同均衡算法的误码率曲线图 这张图是怎么来的呢？运行下面的代码即可。

%%%%%%%%%%%%%%%%%%%%% OFDM仿真 %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% ofdm_fading_sim4.m  %%%%%%%%%
%%%%%%%%%  data:2020年11月17日  author:飞蓬大将军 %%%%%%%%%%


%%%%%程序说明
%%%多径衰落信道下的OFDM传输
%%%调制方式：QPSK
%%%信道编码方式：无
%%%导频方式：梳状类型

%%%sim系列说明
%%%sim2:接收端不做捕获和同步，也不做信道估计，假设知道信道条件，后续版本考虑信道估计
%%%sim3:发送端插入梳状导频、接收端采用LS和MMSE信道估计算法
%%%sim4:尝试解决sim3中遗留问题，并将LS均衡方法、MMSE均衡方法与完美均衡(即接收端知道h，直接fft来做均衡)

%%%%    仿真环境
%软件版本：MATLAB R2019a

%********************** 程序主体 ************%

%%%%%%%%%%%%%%%%%%%%%  参数设置   %%%%%%%%%%%%%%%%%%%
para = 48;   %Number of parallel channel to transmit
fftlen = 64;  %FFT length
noc = 64;    %Number of carrier 
nd = 1;  %Number of information OFDM symbol for one loop
ml = 2;   %Modulation：QPSK
sr = 250000;  %Symbol rate 符号速率
br = sr.*ml;  %Bit rate per carrier
gilen = 16; %length of guard interval 

%%%%%%%%%%比特信噪比设置
%%%先设置大的，程序运行正确后，再设更多的信噪比设置，以画误码率曲线
% ebn0_temp = 80;
ebn0_temp = 1:1:30;
ber_fading_ls_linear = zeros(1,length(ebn0_temp));
ber_fading_ls_spline = zeros(1,length(ebn0_temp));
ber_fading_mmse = zeros(1,length(ebn0_temp));
ber_fading_per = zeros(1,length(ebn0_temp));

%%%%%%%导频信息
Nps = 4;  %导频间隔，起始
B_pilot = 1;  %导频起始位置，第1号子载波
Np = fftlen/Nps; %导频数量

for kkk = 1:length(ebn0_temp)
    
ebn0 = ebn0_temp(kkk);  %Eb/No 


%%%%%%%%%%%%%%%%%%%%%%%%%%  Fading initialization %%%%%%%%%%%%%%

PowerdB=[-1 -8 -17 -21 -25]; % 信道抽头功率特性
Delay=[0 3 5 6 8];          % 信道时延,示例
% Delay=[0 3 5 56 78];          % 信道时延
Power=10.^(PowerdB/10);     % 信道抽头功率特性 '线性'
Ntap=length(PowerdB);       % 信道抽头数
Lch=Delay(end)+1;           % 信道长度


%%%%%%%%%%%%%%%%%%%%%  主循环 %%%%%%%%%%%%%

nloop = 20000; %Number of sumulation loops 

noe2_ls_linear_temp = 0; %Number of error data  of LS_linear
noe2_ls_spline_temp = 0; %Number of error data of  LS_spline
noe2_mmse_temp = 0; %Number of error data of MMSE
noe2_per_temp = 0; %Number of error data of perfect compensation


nod = 0; %Number of transmitted data 

for iii = 1:nloop
    %%%%%%%%%%%%%%%%%  发射机  %%%%%%%%%%%%%%%%%%%
   seldata = rand(1,para*nd*ml)>0.5;  %串行数据

%     seldata = ones(1,para*nd*ml);  %用于调试程序
    paradata = reshape(seldata,para,nd*ml); %串并转换

    [ich,qch] = qpskmod(paradata,para,nd,ml); %调制
    kmod = 1/sqrt(2);
    ich1 = ich.*kmod;
    qch1 = qch.*kmod;
    ch1 = ich1 + qch1*1j;
    
    %%%%%%%%%%%将导频进行插入
%     Xp = 2*(randn(1,Np)>0)-1;    % Pilot sequence generation
%     Xp = 2*randi([0 1],1,Np) - 1;
    Xp = 2*ones(Np,nd) - 1;
    %%%方式一：
    X = zeros(fftlen,nd);
    ip = 0;    
    pilot_loc = [];
    for k=1:fftlen
       if mod(k,Nps)==1
             temp = floor(k/Nps)+1;
             X(k,:) = Xp(temp,:);
%            X(k) = Xp(floor(k/Nps)+1); 
             pilot_loc = [pilot_loc k]; 
             ip = ip+1;
         else
             X(k,:) = ch1(k-ip,:);
       end
    end

      %%%%%方式二：可以先把导频位置和数据位置算出来，比如kk1，kk2，相应放入
    
    %%%%%%%%%%%% IFFT %%%%%%%%%%

    y = ifft(X);
    ich2 = real(y);
    qch2 = imag(y);
    
     %%%观察信号功率
%     spow2 = sum(sum(ich2.^2+ qch2.^2))/nd./para;
    
    
    %%%%%%%% 添加保护间隔 %%%%%%%%
    [ich3,qch3] = giins(ich2,qch2,fftlen,gilen,nd);
    fftlen2 = fftlen + gilen;
    
    
    %****************  Attenuation Calculation **************
%     %%%方式一：
%     spow = sum(ich3.^2+ qch3.^2)/nd./para;
%     attn = 0.5*spow*sr/br*10.^(-ebn0/10);
%     attn = sqrt(attn);
   
    
    %%%方式二：
%     snr = ebn0 + 10*log10(2);
%     attn = sqrt(10.^(-snr/10)*spow/2);
    %以上两种方式表达是一样的
    
    %%%方式三：
    
    spow4 = sum(ich3.^2+ qch3.^2)/nd./fftlen2;
    esn0 = ebn0 + 10*log10(para/fftlen2) + 10*log10(ml);
    attn2 = 0.5*spow4*10.^(-esn0/10);
    attn2 = sqrt(attn2);
    
 
    
%     aaa = 1; %用于调试

    %***************  衰落信道 Fading channel ***************%
    channel = (randn(1,Ntap) + 1j * randn(1,Ntap)).*sqrt(Power/2);
    h = zeros(1,Lch);
    h(Delay+1) = channel;
    y = conv(ich3 + 1j*qch3,h);
    ifade = real(y(:,1:length(ich3)));
    qfade = imag(y(:,1:length(ich3)));
    spow5 = sum(ifade.^2+ qfade.^2)/nd./fftlen2;
    
%     %%%方式四：
%     esn0 = ebn0 + 10*log10(para/fftlen2) + 10*log10(ml);
%     attn2 = 0.5*spow5*10.^(-esn0/10);
%     attn2 = sqrt(attn2);
%     
    %%%%%不管哪种方式，以attn传入
    attn =  attn2;
    
    
    %***********************  接收机 *******************%
    %%%%%%%%%%假设已经完美同步上
    %%%%%%%%%% AWGN addition %%%%%%%%%
    [ich4,qch4] = comb(ifade,qfade,attn);

    %%%%%%%% 去掉保护间隔 %%%%%%%%
    [ich5,qch5] = girem(ich4,qch4,fftlen2,gilen,nd);
    
    
    %%%%%%%%%%%%%%  FFT %%%%%%%%
    rx = ich5 + qch5.*1i;
    ry = fft(rx);
%     ich6 = real(ry);
%     qch6 = imag(ry);

    %%%%%%%  信道估计  %%%%
    for m=1:3
         if m==1 
             H_est_ls_linear = LS_CE(ry.',Xp.',pilot_loc,fftlen,Nps,'linear');
%              method='LS-linear'; % LS estimation with linear interpolation
         elseif m==2
             H_est_ls_spline = LS_CE(ry.',Xp.',pilot_loc,fftlen,Nps,'spline'); 
%              method='LS-spline'; % LS estimation with spline interpolation
         else
%              H_est_mmse = MMSE_CE(ry.',Xp.',pilot_loc,fftlen,Nps,h,SNR); 
              H_est_mmse = MMSE_CE(ry.',Xp.',pilot_loc,fftlen,Nps,h,ebn0); 
%              method='MMSE'; % MMSE estimation
         end
    end % end for count'
    
   
    
    %%%%%%%%  信道均衡 %%%%
    %%%注意A的共轭转置和转置的区别，前者是A'，后者是A.'
    %%%%采用导频估计出的h，进行均衡
    ry_ls_linear_temp = ry./(H_est_ls_linear.');
    ry_ls_spline_temp = ry./(H_est_ls_spline.');
    ry_mmse_temp = ry./(H_est_mmse.');
    
    %%%假设接收端完美知道h，进行均衡
    H = fft([h,zeros(1,fftlen-Lch)].');
    ry_per_temp = ry./H; 
    
    
    
    
    %%%%%%%%%% 去除导频
    ip = 0;
    for k=1:fftlen
       if mod(k,Nps)==1
%              temp = floor(k/Nps)+1;
%              X(k,:) = Xp(temp,:);
% %            X(k) = Xp(floor(k/Nps)+1); 
%              pilot_loc = [pilot_loc k]; 
             ip = ip+1;
       else
             ry_ls_linear(k-ip,:) = ry_ls_linear_temp(k,:);
             ry_ls_spline(k-ip,:) = ry_ls_spline_temp(k,:);
             ry_mmse(k-ip,:) = ry_mmse_temp(k,:);
             ry_per(k-ip,:) = ry_per_temp(k,:);
       end
    end
    
    
    %%%%%%%%%%%%%% demoluation %%%%%%%%%%%%%
    demodata_ls_linear = qpskdemod(real(ry_ls_linear)./kmod,imag(ry_ls_linear)./kmod,para,nd,ml);
    demodata_ls_spline= qpskdemod(real(ry_ls_spline)./kmod,imag(ry_ls_spline)./kmod,para,nd,ml);
    demodata_mmse = qpskdemod(real(ry_mmse)./kmod,imag(ry_mmse)./kmod,para,nd,ml);
    demodata_per = qpskdemod(real(ry_per)./kmod,imag(ry_per)./kmod,para,nd,ml);
    
    %%%%%%%%%%%%%% 并串转换 %%%%%%%%%
    demodata1_ls_linear = reshape(demodata_ls_linear,1,para*nd*ml);
    demodata1_ls_spline = reshape(demodata_ls_spline,1,para*nd*ml);
    demodata1_mmse  = reshape(demodata_mmse,1,para*nd*ml);
    demodata1_per = reshape(demodata_per,1,para*nd*ml);
    
    
    %%%%%%%%%%%%%%% Bit Error Rate %%%%%%%%%%%
    noe2_ls_linear = sum(abs(demodata1_ls_linear-seldata));
    noe2_ls_spline = sum(abs(demodata1_ls_spline-seldata));
    noe2_mmse = sum(abs(demodata1_mmse-seldata));
    noe2_per = sum(abs(demodata1_per-seldata));
    nod2 = length(seldata);
    
    %%%%cumulative the number of error and data in noe and nod
    noe2_ls_linear_temp = noe2_ls_linear_temp + noe2_ls_linear; %Number of error data  of LS_linear
    noe2_ls_spline_temp = noe2_ls_spline_temp + noe2_ls_spline ; %Number of error data of  LS_spline
    noe2_mmse_temp = noe2_mmse_temp + noe2_mmse; %Number of error data of MMSE
    noe2_per_temp = noe2_per_temp + noe2_per;
    nod = nod + nod2;
        
    
end
%**************   output result *************
ber_ls_linear_temp = noe2_ls_linear_temp/nod;
ber_ls_spline_temp = noe2_ls_spline_temp/nod;
ber_mmse_temp = noe2_mmse_temp/nod;
ber_per_temp = noe2_per_temp/nod;


ber_fading_ls_linear(1,kkk) = ber_ls_linear_temp ;
ber_fading_ls_spline(1,kkk) = ber_ls_spline_temp;
ber_fading_mmse(1,kkk) = ber_mmse_temp;
ber_fading_per(1,kkk) = ber_per_temp;


end

%************************* 画误码率曲线进行对比 ******************%
ebn0_temp = 1:1:30;
rayleign_one_path_theory = ber_temp(ebn0_temp); 
semilogy(ebn0_temp,rayleign_one_path_theory,'-*',ebn0_temp,ber_fading_ls_linear,'-^',ebn0_temp,ber_fading_ls_spline,'->',ebn0_temp,ber_fading_mmse,'-<',ebn0_temp,ber_fading_per,'-+');
xlabel('比特信噪比');
ylabel('误码率');
title('多径衰落信道下误码率仿真曲线');
legend('理论曲线','lslinear实验曲线','lsspline实验曲线','mmse实验曲线','完美均衡实验曲线');
grid on;
%***********   多径瑞利衰落下的OFDM采用BPSK调制的误码率值 **********%

%%%%%%%%%%%%%%%%%%%%%%%      理论值          **************%
%%%%%%%%%%%%%     EbN0(dB)      误码率        
%%%%%%%%%%%%%       3        0.125000000000000
%%%%%%%%%%%%%       4        0.100000000000000
%%%%%%%%%%%%%       5        0.0833333333333333
%%%%%%%%%%%%%       6        0.0714285714285715
%%%%%%%%%%%%%       7        0.0625000000000000
%%%%%%%%%%%%%       8        0.0555555555555556
%%%%%%%%%%%%%       9        0.0500000000000000
%%%%%%%%%%%%%      10        0.0454545454545455

%%%%%%%%%%%%%%%%%   结论    %%%%%%%%%%%%%%%%%%%%
%%%实验记录 2020年11月14日
%含已知导频的OFDM经过多径衰落信道的误码率仿真
%OFDM中引入保护间隔，是一种冗余信息，因此相比于理论误码率曲线有10*log（160/128）=0.969dB的损失
%OFDM信道估计引入了发送和接收端都已知的导频信号，这也是冗余，所以设计SNR与Eb/N0的换算
%本次实验还缺少相应的理论分析

运行上面的代码，完整的代码已经上传github，地址：123kevin456/OFDM-

有以下几点需要注意，当然也是我自己的思考，不一定就完全正确，欢迎你提问、交流与讨论：

（1）在上一讲《OFDM 信道估计与均衡（3）》中讲到过AWGN信道和衰落信道h的不同，对于高斯白噪声信道来说，在发送端和接收端信号功率是没有变化的。但是对于衰落信道来说，衰落信道会使得发送功率和到达接收端功率不一样。

因此，我在上面的代码中分别提到了四种方式，目的是计算并观察信号在不同位置的功率，比如经过IFFT前后，经过衰落信道前后等等。

实际系统中，接收端的白噪声功率是一个固定值。

比如接收端的噪声功率谱密度是0.5w/Hz（也就是0.5J），发送端一个1J的比特发过来。如果是AWGN，那么到达接收端这个比特的能量仍然是1J，此时接收端比特信噪比是3dB。

但是若经过h=0.5的单径衰落信道过来接收端，那么到达接收端这个比特的能量就是0.25J了，此时比特信噪比是-3dB。

所以上面代码中h的设置，其实会影响到图2曲线的左移或者右移。至于是不是正确，仿真是不是应该这样仿，我不太确定，哈哈哈哈哈哈。

（2）仿真中我们都先习惯于设置Eb/N0的值，然后去推Es/N0的值，这个怎么推，我在《第二章：线性分组码》和《第1章：BPSK调制解调器仿真》讲过，即需要考虑到调制编码的影响。

有了采样率和系统带宽后，便可以继续从Es/N0往SNR推。

但在OFDM中有了CP和导频后，这些数据会影响到Eb/N0与Es/N0换算中的k值吗？

我觉得会的，因为CP和导频的存在，本质上都是为了通信系统的可靠性。

CP和导频均消耗了系统的功率资源，但都没有为系统发送更多的数据做贡献，所以将他们当做冗余来进行处理。这也是为什么多次尝试在上面的代码总修改功率的计算位置。

无论是《MIMO-OFDM无线通信技术及MATLAB实现》还是《Simulation and Software Radio for Mobile Communications》，均没有我上面疑问的答案。

在仿真中是否需要我这么考虑，我也不确定，造成的结果也是实验曲线会相比理论有一定的平移。

至此，OFDM信道估计与均衡部分就全部讲完。下一步可能会更新OFDM同步技术。

四、总结

自本科听说过OFDM、信道估计、均衡等概念，到前段时间开始思考如何用MATLAB来仿真OFDM系统呢？相隔了三年。

看到了知乎上@子木的《OFDM完整仿真过程及解释（MATLAB）》以及《给“小白”图示讲解OFDM的原理》，他们的文章写得非常好，帮助我开始理解OFDM相关的知识。

然后我结合《MIMO-OFDM无线通信技术及MATLAB实现》和《Simulation and Software Radio for Mobile Communications》书中的代码，开始慢慢摸索写出仿真。但是在调试代码的过程中，常常出现问题。偶与师兄讨论，有所收获，记录在此，希望对刚接触OFDM仿真的同学有所帮助。

欢迎你双击屏幕、点赞、收藏、转发和分享，关注我的知乎号、CSDN号，也欢迎读者朋友就相关技术问题与我交流，一起学习，共同进步。请你也别忘了把这篇文章分享给你身边正在学习通信专业的同学们，也许能够帮到Ta。

这是《陈老湿·通信MATLAB》仿真的第7章，期待下次更新见！