OFDM同步算法之Park算法

park算法代码

训练序列结构 T=[(C) (D) (C^{*}) (D^{*})],其中C表示由长度为N/4的复伪随机序列PN，ifft变换得到的符号序列
(C(n) = D(N/4-n))

---

原文解释：The training symbol is produced by transmitting
a real-valued PN sequence on the even frequencies, while zeros
are used on the odd frequencies. This means that one of the
points of a BPSK constellation is transmitted at each even fre-quency.

(igstar)park：为了进一步解决minn算法存在的不足，park等人在分析了schmidl算法和minn算法定时效果不佳的原因后，重新设计了新的前导训练序列的结构，并给出了新的定时同步度量函数，该算法的定时度量函数曲线出现了一个更为尖锐的自相关峰，很明显该算法消除了schmidl算法中由于循环前缀的存在而导致的平顶效应，同时得到了比minn算法更为尖锐的自相关峰，提高了定时的精度和确定性，但是在噪声干扰较大的情况下，该算法还是会出现较大的定时同步估计误差，其同步估计的稳定性依然较差。

参考文献
Park B,Choen H , KO E ,et al.A novel timing estimation method for OFDM systems[J].IEEE Commun.Leet.2003,7(5):53-55.

[M(d)=frac{left | P(d) ight |}{R^{2}(d)}^{2} ]

[P(d)=sum_{m=0}^{N/2 -1}r(d-m) r(d+m) ]

[R(d)=sum_{m=0}^{N/2-1}left | r(d+m) ight |^{2} ]

实际在算法实现上
(P(d)=sum_{m=0}^{N/2-1}r(d-1-m) r(d+m))
这是因为序列个数通常是偶数而非奇数，不会出现
(r(d)r(d))的情况。
所求得的d对应的是训练序列（不包含循环前缀）的中间位置。

仿真验证如果发送的数据是随机的[1+1i,-1+1i,-1-1i,1-1i],且训练队列由PN序列（用随机序列代替）通过IFFT得到时，结果与原论文结果最相近。而且如果此时PN序列的最大值为7时更容易看到跟原论文作者一样的结果，具有两个副峰

clear all; 
clc; 
%参数定义 
N=256;       %FFT/IFFT 变换的点数或者子载波个数（Nu=N） 
Ng=N/8;      %循环前缀的长度 (保护间隔的长度) 
Ns=Ng+N;     %包括循环前缀的符号长度 
SNR=25;
%************利用查表法生成复随机序列********************** 
QAMTable=[7+7i,-7+7i,-7-7i,7-7i]; 
buf=QAMTable(randi([0,3],N/2,1)+1); %加1是为了下标可能是0不合法
 
%产生train 
pn=rand(1,N/2)>0.5; 
pn=reshape(pn,N/4,2); 
[ich,qch]=qpskmod(pn,N/4,1,2); 
kmod=sqrt(2); 
x=ich*kmod+qch*kmod*i; 
y=ifft(x); 
y=reshape(y,N/4,1); 
train=[y;y(N/4:-1:1,1);conj(y);conj(y(N/4:-1:1,1))];
 
%*****************添加一个空符号以及一个后缀符号************* 
src = QAMTable(randi([0,3],N,1)+1).'; 
sym = ifft(src); 
sig =[zeros(N,1) train sym]; 
%sig =[sym train sym];
 
%**********************添加循环前缀************************* 
tx =[sig(N - Ng +1:N,:);sig];
%tx = [sig(1,N-Ng+1:N) sig];
 
%***********************经过信道*************************** 
recv = reshape(tx,1,size(tx,1)*size(tx,2)); %size的1表示行，2表示列，从%前向后数，超过了为1
%recv = tx;
%recv1 = awgn(recv,1,'measured'); 
%recv2 = awgn(recv,5,'measured'); 
%recv3 = awgn(recv,10,'measured');
%recv = awgn(recv,SNR);
%*****************计算符号定时***************************** 
P=zeros(1,2*Ns); 
R=zeros(1,2*Ns); 

for d = Ns/2+1:1:2*Ns 
    for m=0:N/2  
        P(d-Ns/2) = P(d-Ns/2) + (recv(d+m))*recv(d-1-m);  
        R(d-Ns/2) = R(d-Ns/2) + power(abs(recv(d+m)),2); 
    end 
end 
% for d = Ns/2+1:1:2*Ns 
%     for m=0:1:(N/2-1)  
%         P(d-Ns/2) = P(d-Ns/2) + recv(d-m)*recv(d+m);  
%         R(d-Ns/2) = R(d-Ns/2) + power(abs(recv(d+m)),2); 
%     end 
% end 
M=power(abs(P),2)./power(abs(R),2); 
[a b]=max(M);
b+Ns/2
%**********************绘图****************************** 
figure('Color','w'); 
d=1:1:400; 
figure(1);
plot(d,M(d+N/2)); 
grid on; 
axis([0,400,0,1.1]); 
title('park algorithm'); 
xlabel('Time (sample)'); 
ylabel('Timing Metric'); 
%legend('no noise','SNR=1dB','SNR=5dB','SNR=10dB');)); 
hold on;