初认OFDM（四）：CFO的影响和估计

一.什么是CFO

这个问题其实在上一张博客内已经解释过了，此处直接复制过来

STO是符号定时偏差，它指的是OFDM中FFT窗未能完全与符号未能精确对齐而导致的误差，STO和CFO是OFDM同步技术的两个方向，
$\begin{aligned} y_l[n]&=IDFT\{Y_l[k]\}=IDFT\{H_l[k]X_l[k]+Z_l[k]\} \\ &=\frac{1}{N}\sum_{k=0}^{N-1}H_l[k]X_l[k]e^{j2\pi (k+\epsilon)(n+\delta)/N}+z_l[n] \end{aligned}$
我是这样理解这个式子的，还记得上面的式子里k和n是如何来的吗，k来自于子载波的频率采样， $f_k=k/T_{sym}$ ,因此CFO就加在了k上，而n来自 $t=lT_{sym}+nT_s$ ,因此它对应STO

二. CFO的影响

CFO主要来自与两个方面

本地振荡器的频率不统一而带来的无法相干，这也是波束赋形最主要想解决的问题
多普勒频移

和STO一样，CFO对接受信号的影响如图

但请注意，陈老湿的博客已经指出，这个表格是错误的，应该是 $x[n]e^{j2\pi n \epsilon/N}$ ,这表示，时域的OFDM符号是有相位偏移的，而且随着n增大他也在增大，注意这个n并不单单指一个符号，而是整个OFDM序列

三. IFO和CFO

1.IFO

IFO相差了整数倍，不会影响ICI，但其实也难以估计，因为你可能压根不知道它的存在

2.FFO的影响

四 CFO估计

CFO的估计思路是类似的

使用CP，CP与符号末尾应该是一致的，那么每个 $x [n]$ 相差的相位可以计算
使用训练符号内部重复，那么每个 $x [n]$ 相差的相位可以计算
CFO的问题当然要放到频域来处理，发送两个一样的训练符号，那么他们在频域提取公因式应该相差一个相位（Moose）
Classen方法目前我还没弄懂，看明白后会更新

五. OFDM估matlab代码

%完成时域基于CP的方法和频域的Moose/Classen方法，用于后续CFO补偿
clear, clf
CFO = 0.15;                     % CFO大小
% CFO = 0;
Nfft=128;                       % FFT采样数  
Nbps=2;                         % QPSK或QAM
M=2^Nbps;                       % 每个符号代表几比特
Es=1;                           % 能量
A=sqrt(3/2/(M-1)*Es);           % QAM归一化
N=Nfft;                         % 发送的符号长度，为了方便，和Nfft保持一致
Ng=Nfft/4;                      % GI长度
Nofdm=Nfft+Ng;                  % 一个OFDM符号的长度
Nsym=3;                         % 一共发送了3个OFDM符号，前两个是导频，最后一个是真正发送的数据符号
% h=complex(randn,randn)/sqrt(2);
% %h=[1 zeros(1,5)]; 
% channel(h,0);  
%Transmit signal
x=[];
for m=1:Nsym                    % 前两个是导频，最后一个是真正发送的数据符号
   msgint=randi([0 M-1],1,N);   % 生成要发送的符号

   if m<=2                      % 
       Xp = add_pilot(zeros(1,Nfft),Nfft,4);    % 生成导频
       Xf=Xp; % add_pilot
   else  %Xf= QAM(msgint((i-1)*N+1:i*N),Nbps);  % constellation mapping. average power=1        
       Xf = A.*qammod(msgint,M,'UnitAveragePower',true);
   end                                      
   xt = ifft(Xf,Nfft);          % ifft  
   x_sym = add_CP(xt,Ng);       % 加CP
   x= [x x_sym];                % 加CP
end    

%channel 可添加所需信道
y=x; % No channel effect

%Signal power calculation
sig_pow= y*y'/length(y); % Signal power calculation

%%%%
SNRdBs= 0:3:30;  
% SNRdBs= 100; 设100是为调试程序  
MaxIter = 100;  
for i=1:length(SNRdBs)
   SNRdB = SNRdBs(i);
   MSE_CFO_CP = 0; 
   MSE_CFO_Moose = 0; 
   MSE_CFO_Classen = 0;
   rand('seed',1);              % 设置种子来保证每次仿真结果一致
   randn('seed',1);
   y_CFO= add_CFO(y,CFO,Nfft);  % 增加CFO，此处是在时域添加的，因此是×相位
   % 多次迭代取平均
   for iter=1:MaxIter
      %y_aw=add_AWGN(y_CFO,sig_pow,SNRdB,'SNR',Nbps);  % AWGN added, signal power=1
      y_aw = awgn(y_CFO,SNRdB,'measured');              % 增加高斯白噪声
      
      Est_CFO_CP = CFO_CP(y_aw,Nfft,Ng); % CP-based     % 根据CP测算CFO
      MSE_CFO_CP = MSE_CFO_CP + (Est_CFO_CP-CFO)^2;     % 平方累计
      
      Est_CFO_Moose = CFO_Moose(y_aw,Nfft);             % Moose估计
      MSE_CFO_Moose = MSE_CFO_Moose + (Est_CFO_Moose-CFO)^2;% 平方累计
      
      Est_CFO_Classen = CFO_Classen(y_aw,Nfft,Ng,Xp); % Classen (Pilot-based)
      MSE_CFO_Classen = MSE_CFO_Classen + (Est_CFO_Classen-CFO)^2;
      
   end % the end of for (iter) loop
   MSE_CP(i) = MSE_CFO_CP/MaxIter; 
   MSE_Moose(i) = MSE_CFO_Moose/MaxIter;  
   MSE_Classen(i) = MSE_CFO_Classen/MaxIter;
end%ebn0 end    
semilogy(SNRdBs, MSE_CP,'-+');
grid on, hold on
semilogy(SNRdBs, MSE_Moose,'-x'); semilogy(SNRdBs, MSE_Classen,'-*');
xlabel('SNR[dB]'), ylabel('MSE'); title('CFO Estimation'); %axis([0 30 10e-8 10e-2])
% str=sprintf('CFO = %1.2f',CFO);
legend('CP-based technique','Moose (Preamble-based)','Classen (Pilot-based)');
% legend(str);

function xp=add_pilot(x,Nfft,Nps)
% CAZAC (Constant Amplitude Zero AutoCorrelation) sequence --> pilot
% Nps : Pilot spacing
% 增加导频，增加的是大名鼎鼎的Zadoff –Chu序列，他是CAZAC序列的一种
% 这个函数输入是一个zeros矩阵
% 这里的Nps不要和主函数里QAM的Nps混为一谈
% 请注意这个add_pilot是在频域生成的，因此Nps就是子载波的间隔，这就是梳状导频
% 这里的Nps是为了让导频只占用部分的子载波，平衡频谱的利用率和估计的准确性
if nargin<3
    Nps=4; 
end
Np=Nfft/Nps;  
xp=x; % Number of pilots and an OFDM signal including pilot signal
for k=1:Np
   xp((k-1)*Nps+1)= exp(1i*pi*(k-1)^2/Np);  % Pilot boosting with an even Np
   %xp((k-1)*Nps+1)= exp(j*pi*(k-1)*k/Np);  % Pilot boosting with an odd Np
end

function CFO_est=CFO_CP(y,Nfft,Ng)
% Time-domain CFO estimation based on CP (Cyclic Prefix)
% 用CP在时域进行估计，因为振荡器或者是多普勒频移，通过y结算Y产生了误差
% 而我们可以把它看作是由于y没有接受到位产生的
% 这就是y产生了相位偏差的来源
% 这个相位和\epsilon还有n是成正比的
% 因此我们可以用CP和OFDM尾部进行一个估计
nn=1:Ng; 
CFO_est = angle(y(nn+Nfft)*y(nn)')/(2*pi);

function CFO_est=CFO_Moose(y,Nfft)
% Frequency-domain CFO estimation using Moose method based on two consecutive identical OFDM symbols 
% 大名鼎鼎的Moose估计
% 这是在频域完成的，他发射了两个相同的导频信号
% 既然相同而之间有相差了N，那么第一个导频的时域的x[n]和第二个导频的x[n]就都相差了2\pi\epsilon
% 做一个DFT后，提取相同的系数，当然频域也相差了这么多，这样就可以估计了
for i=0:1   
   Y(i+1,:)= fft(y(Nfft*i+1:Nfft*(i+1)),Nfft);
end
CFO_est = angle(Y(2,:)*Y(1,:)')/(2*pi);

function CFO_est=CFO_Classen(yp,Nfft,Ng,Nps)
% Frequency-domain CFO estimation using Classen method based on pilot tones
% 频域估计，Classen估计法
% 计算的是两个连续的OFDM符号中的导频信号之间的相位差
% 总之生成一个导频
if length(Nps)==1
    Xp=add_pilot(zeros(1,Nfft),Nfft,Nps); 
else
    Xp=Nps; % If Nps is given as an array, it must be a pilot sequence Xp
end
Nofdm=Nfft+Ng; 
% 找出导频的子载波
kk=find(Xp~=0); 
Xp=Xp(kk); % Extract pilot tones
for i=1:2 
   yp_without_CP = remove_CP(yp(1+Nofdm*(i-1):Nofdm*i),Ng);
   Yp(i,:) = fft(yp_without_CP,Nfft);
end
CFO_est = angle(Yp(2,kk).*Xp*(Yp(1,kk).*Xp)')/(2*pi)*Nfft/Nofdm; % Eq.(5.31)