初识OFDM（三）：OFDM同步技术之STO

一.什么是STO

首先要弄明白，STO是发生在FFT中的，而不是IFFT，也就是说，他在接受时出现了问题

STO是符号定时偏差，它指的是OFDM中FFT窗未能完全与符号未能精确对齐而导致的误差，STO和CFO是OFDM同步技术的两个方向，CFO我们会在下一节介绍，不过他们可以用同一个式子来表示
y l [ n ] = I D F T { Y l [ k ] } = I D F T { H l [ k ] X l [ k ] + Z l [ k ] } = 1 N ∑ k = 0 N − 1 H l [ k ] X l [ k ] e j 2 π ( k + ϵ ) ( n + δ ) / N + z l [ n ] \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} yl​[n]​=IDFT{Yl​[k]}=IDFT{Hl​[k]Xl​[k]+Zl​[k]}=N1​k=0∑N−1​Hl​[k]Xl​[k]ej2π(k+ϵ)(n+δ)/N+zl​[n]​
我是这样理解这个式子的，还记得上面的式子里k和n是如何来的吗，k来自于子载波的频率采样，$f_k=k/T_{sym}$,因此CFO就加在了k上，而n来自$t=l{T}_{sym}+n{T}_s$,因此它对应STO

二. STO的影响

时域上的偏差导致了频域上相位的偏差。

而STO的影响可以根据下图分为四类

1. case1

精确定时，没有干扰

2.case2

这种情况是提前了，但是没有被上一个符号的时延扩展影响，这样是不会有ISI的,注意这里的$\delta$实际上是负的，此外书上应该写错了，应该是对$y_l$做FFT

3 case3

4 case4

落后于正常时间，这时候，由于CP的存在，FFT窗中实际上混入了下一个OFDM符号

理解一下，CP正常的位置是$x_{l+1}[n+\delta -N_G]$,再偏移$$\delta$$

5 matlab实现

这里的实现有两点需要指出

1. y这里指的是一个由多个OFDM符号组成的OFD吗帧，该函数的作用仅仅是平移，而不是复制CP
2. 关于nSTO的正负，请看我在后续代码中的注释

三. STO估计技术

1时域CP估计

CP估计思路如下图

既然现在接受的不对，那我微调开始处，找到CP

这里说明一下，本质上都是比较W1窗和W2窗的相似性，比较使用了两种方法

• 比较差，差最小时最大，并有一些修正方法
• 比较相关性，相关性最高时最大，并有一些修正方法

注意，这里的示意图画的是有一些问题的，窗口是从$$\delta$$处开始的，结束于CP结尾处，而不是和CP对其的

2 训练符号的STO估计

训练符号实际上就是前导，发送已知的符号，在数据符号中间隔插入

书上说的不是很清楚，我的理解是：可以看到，训练符号中的CP并不是尾部的重复而是不一样的数据，W1和W2同时滑动时，和两个C分别重合时获得最相似，进而估计

3 频域STO估计

四.代码详解

需要指出，这里的代码仍然来自于

%基于CP，采用最大相关和最小距离算法完成STO的估计

clear, figure(1), clf, figure(2), clf
nSTOs = [-3 -3 2 2];                %对应STO的采样数
CFOs = [0 0.5 0 0.5];               %CFO向量
SNRdB = 30;                         %SNR
MaxIter = 10;                       %迭代次数
%CFOs = [0 0 0 0];
Nfft = 128;                         %FFT大小
Ng = Nfft/4;                        % GI长度
Nofdm = Nfft+Ng;                    % OFDM符号长度
Nbps = 2;                           % 2/4 对应 QPSK/16QAM 
M = 2^Nbps;                         % 符号对应的可能性数量
Es = 1; 
A = sqrt(3/2/(M-1)*Es);             % QAM归一化因子
N = Nfft; 
com_delay = Nofdm/2;                % 公共时延
Nsym = 100;                         % 一共有一百个OFDM符号
rand('seed',1);                     % 设置种子
randn('seed',1);
for i=1:length(nSTOs)               % 对于每一个不同的\delta STO
   nSTO= nSTOs(i);  
   CFO= CFOs(i);
   x = [];                          % 初始化信号块，最后x就是发送出来的OFDM模块
   for m=1:Nsym % random bits generates 
      msgint=randi([0 M-1],1,N);    % 生成传输符号,1*128个0-3的数
      Xf = A.*qammod(msgint,M,'UnitAveragePower',true);% 调制成复数
      %***********************缺少了倒置的过程***********************%
      xt = ifft(Xf,Nfft);           % 发送
      x_sym = add_CP(xt,Ng);        % 加CP
      x = [x x_sym];
   end
   %*********************** 信道 ************************%
   %%%%%在这里根据需求添加信道，先假设是没有信道
   y = x;  % 没有信道影响
   
   sig_pow = y*y'/length(y);        % 计算能量，sig_pow= mean(mean(y.*conj(y),2))
   
   % 频率偏移 +  符号定时偏移 
   y_CFO= add_CFO(y,CFO,Nfft);      % 加CFO
   y_CFO_STO= add_STO(y_CFO,-nSTO); % 加STO，这是加在整个信号上的，因此头尾补零就行了，但是取-是因为？在研究
   v_ML=zeros(1,Ng);                % 初始化
   v_Cl=zeros(1,Ng);
   Mag_cor= 0;                      % arg的结果
   Mag_dif= 0;
   %%Add additive white gaussian noise
   for iter=1:MaxIter
      % 加噪声
      y_aw = awgn(y_CFO_STO,SNRdB,'measured');
      
      %%%%%%%Symbol Timing Acqusition
      [STO_cor,mag_cor]= STO_by_correlation(y_aw,Nfft,Ng,com_delay); %书中自带
      [STO_cor_temp,mag_cor_temp]= STO_by_correlation_sim1(y_aw,Nfft,Ng,com_delay); %我自己编写
      %%%%%经验证，以上两者函数结果一致
      
      
      [STO_dif,mag_dif] = STO_by_difference(y_aw,Nfft,Ng,com_delay); %书中自带
      [STO_dif_temp,mag_dif_temp] = STO_by_difference_sim1(y_aw,Nfft,Ng,com_delay); %我自己编写
      %%%%%经验证，以上两者函数结果一致
      % 计数
      % 这里取反了，返回了符合"左加右减"的直觉的STO
      v_ML(-STO_cor+Ng/2)= v_ML(-STO_cor+Ng/2)+1;
      v_Cl(-STO_dif+Ng/2)= v_Cl(-STO_dif+Ng/2)+1;
      Mag_cor= Mag_cor + mag_cor;
      Mag_dif= Mag_dif + mag_dif;
   end % End of for loop of iter
   %%%%%%% Probability
   v_ML_v_Cl = [v_ML; v_Cl]*(100/MaxIter);      % 取百分数
   figure(1+i-1); 
   set(gca,'fontsize',9);
%    subplot(220+i)
    
   bar(-Ng/2+1:Ng/2,v_ML_v_Cl');
   hold on, grid on
   str = sprintf('nSTO Estimation: nSTO=%d, CQFO=%1.2f, SNR=%3d[dB]',nSTO,CFO,SNRdB);           
   title(str); 
   xlabel('Sample'), ylabel('Probability');
   legend('ML','Classen'); 
   axis([-Ng/2-1 Ng/2+1 0 100])
   %%%%%%% Time metric
   Mag_cor = Mag_cor/MaxIter; 
   [Mag_cor_max,ind_max] = max(Mag_cor);
   nc= ind_max-1-com_delay;
   Mag_dif = Mag_dif/MaxIter; 
   [Mag_dif_min,ind_min] = min(Mag_dif);
   nd= ind_min-1-com_delay
   nn= -Nofdm/2 + [0:length(Mag_cor)-1]; 
   nt= nSTO;
%    figure(2);
%    subplot(220+i);
   figure(5+i-1); 
   plot(nn,Mag_cor,nn,1.5*Mag_dif,'r:','markersize',1);
   hold on
     
   stem(nc,Mag_cor_max,'b','markersize',5);
   stem(nSTO,Mag_cor(nSTO+com_delay+1),'k.','markersize',5); % Estimated/True Maximum value
   str1 = sprintf('STO Estimation - ML(b-)/Classen(r:) for nSTO=%d, CFO=%1.2f',nSTO,CFO); %,SNRdB);
   title(str1); 
   xlabel('Sample'), ylabel('Magnitude'); 
 
   %stem(n1,Mag_dif_min,'r','markersize',5)
   stem(nd,Mag_dif(nd+com_delay+1),'r','markersize',5);
   stem(nSTO,Mag_dif(nSTO+com_delay+1),'k.','markersize',5); % Estimated/True Minimum value
   set(gca,'fontsize',9, 'XLim',[-32 32], 'XTick',[-10 -3 0 2 10]); %, xlim([-50 50]),
   legend('基于相关的','基于差值最小的'); 
end % End of for loop of i

function [STO_est,Mag] = STO_by_difference(y,Nfft,Ng,com_delay)
%通过最小化 CP和 OFDM 符号后部的差值，实现STO的估计
%输入:
%y:包括 CP的OFDM接收信号
%Ng: CP/GI内的采样数
%com_delay:公共时延
%输出:
%STO est :STO估计
%Mag:相关函数的时变轨迹
N_ofdm = Nfft+Ng;
minimum = 100;                          % 初始化最小值
STO_est = 0;                            % 初始化估计结果
if nargin<4                             
    com_delay = N_ofdm/2;
end
% 这里仔细考虑一下，为什么STO是这样计算的
% 如果STO为0，即没有偏移，那(n+com_delay)=等于N_ofdm的整数倍
% 现在如果STO是一个大于零的值如3，就是说接收到的y整体滞后了，CP头本应该是N_ofdm，现在变成了N_ofdm-3
% 那最终(n+com_delay)=N_ofdm-3，最终得到的STO就是3
% 这也就是之前取反的原因
% 这其实是不符合直觉"左加右减"的，FFT窗滞后了，\delta是一个正值
for n = 1:N_ofdm
    nn = n+com_delay+[0:Ng-1];          % 这就是窗口了，公共时延就是窗开始的位置，这里是符号的中间
    tmp0 = abs(y(nn))-abs(y(nn+Nfft));  % 计算差值
    Mag(n) = tmp0*tmp0';                % 由式(5.11)给出的差值的平方
    if Mag(n)< minimum                  
        minimum = Mag(n);               % 找到最小值
        STO_est = N_ofdm-com_delay -(n-1);  % 最后的估计结果，符号长度-公共实验-当前位置
    end
end

function [STO_est, Mag] = STO_by_correlation(y,Nfft,Ng,com_delay)%通过最大化CP和 OFDM符号后部的相关函数,得到STO的估计%输入:
%   y:包括CP的 OFDM接收信号
%Ng: GICP内的采样数
%com_delay :公共时延
%输出;
%STO_est :STO估计
% Mag:相关函数的时变轨迹
% 该函数和difference是类似的，不多赘述
N_ofdm = Nfft+Ng; % OFDM符号长度
if nargin<4
    com_delay = N_ofdm/2;
end
nn =0:Ng-1;
yy = y(nn+com_delay)*y(nn+com_delay+Nfft)';
%相关函数
%这里的n1才对,difference中的是有问题的
maximum = abs(yy);
for n=1:N_ofdm
    n1 = n-1;
    yy1 = y(n1+com_delay)*y(n1+com_delay+Nfft)';
    yy2 = y(n1+com_delay+Ng)*y(n1+com_delay+Nfft+Ng)';
    yy = yy-yy1+yy2;Mag(n)=abs(yy); %式（5.13）
    if(Mag(n) > maximum)
        maximum = Mag(n);
        STO_est = N_ofdm-com_delay-n1;
    end
end

5 reference

[1] 通信原理 第2版 [李晓峰 编著] 2014年版

[2] MIMO-OFDM 无线通信技术及MATLAB 实现

[3] https://zhuanlan.zhihu.com/p/389034302