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;

OFDM同步算法之Park算法