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.
---

\(\bigstar\)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) \right |}{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) \right |^{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;