%********************schmidl algorithm******************* 
%ofdm符號(hào)同步算法-schmidl&cox算法的matlab實(shí)現(xiàn)%
%********************************************************* 
close all; 
clear all; 
clc; 
%參數(shù)定義 
N=256;       %FFT/IFFT 變換的點(diǎn)數(shù)或者子載波個(gè)數(shù)（Nu<=N） 
Ng=N/8;      %循環(huán)前綴的長(zhǎng)度 (保護(hù)間隔的長(zhǎng)度) 
Ns=Ng+N;     %包括循環(huán)前綴的符號(hào)長(zhǎng)度 
deltaf=0.45;  %頻偏  
%************利用查表法生成復(fù)隨機(jī)序列********************** 
QAMTable=[7+7i,-7+7i,-7-7i,7-7i]; 
buf=QAMTable(randint(N/2,1,4)+1); 
 
%*************在奇數(shù)子載波的位置插入零********************* 
x=zeros(N,1); 
index = 1; 
for n=1:2:N 
     x(n)=buf(index); 
     index=index+1; 
end; 
 
%**************利用IFFT變換生成Schmidl訓(xùn)練符號(hào)*************** 
sch = ifft(x);   %[A A]的形式 
 
%*****************添加一個(gè)空符號(hào)以及一個(gè)后綴符號(hào)************* 
src = QAMTable(randint(N,1,4)+1).'; 
sym = ifft(src); 
sig =[zeros(N,1) sch sym]; 
 
%**********************添加循環(huán)前綴************************* 
tx =[sig(N - Ng +1:N,:);sig]; 
 
%***********************經(jīng)過信道*************************** 
recv = reshape(tx,1,size(tx,1)*size(tx,2)); 
recv1 = awgn(recv,1,'measured'); 
recv2 = awgn(recv,5,'measured'); 
recv3 = awgn(recv,10,'measured'); 
recv4 = awgn(recv,10,'measured');
for k=1:length(recv4); 
     recv4(k)=recv4(k)*exp(i*2*pi*deltaf*k/N); 
end 

%*****************計(jì)算符號(hào)定時(shí)***************************** 
P=zeros(1,2*Ns); 
R=zeros(1,2*Ns); 
P1=zeros(1,2*Ns); 
R1=zeros(1,2*Ns); 
P2=zeros(1,2*Ns); 
R2=zeros(1,2*Ns); 
P3=zeros(1,2*Ns); 
R3=zeros(1,2*Ns); 
P4=zeros(1,2*Ns); 
R4=zeros(1,2*Ns); 
for d = Ns/2+1:1:2*Ns 
    for m=0:1:N/2-1  
        P(d-Ns/2) = P(d-Ns/2) + conj(recv(d+m))*recv(d+N/2+m);  
        R(d-Ns/2) = R(d-Ns/2) + power(abs(recv(d+N/2+m)),2); 
        P1(d-Ns/2) = P1(d-Ns/2) + conj(recv1(d+m))*recv1(d+N/2+m);  
        R1(d-Ns/2) = R1(d-Ns/2) + power(abs(recv1(d+N/2+m)),2); 
        P2(d-Ns/2) = P2(d-Ns/2) + conj(recv2(d+m))*recv2(d+N/2+m);  
        R2(d-Ns/2) = R2(d-Ns/2) + power(abs(recv2(d+N/2+m)),2); 
        P3(d-Ns/2) = P3(d-Ns/2) + conj(recv3(d+m))*recv3(d+N/2+m);  
        R3(d-Ns/2) = R3(d-Ns/2) + power(abs(recv3(d+N/2+m)),2); 
        P4(d-Ns/2) = P4(d-Ns/2) + conj(recv4(d+m))*recv4(d+N/2+m);  
        R4(d-Ns/2) = R4(d-Ns/2) + power(abs(recv4(d+N/2+m)),2); 
    end 
end 
M=power(abs(P),2)./power(abs(R),2); 
M1=power(abs(P1),2)./power(abs(R1),2); 
M2=power(abs(P2),2)./power(abs(R2),2); 
M3=power(abs(P3),2)./power(abs(R3),2); 
M4=power(abs(P4),2)./power(abs(R4),2); 
a=max(P4); % 求最大值
a1=find(P4==a);% 求最大值所在位置 
f=angle(P4(a1))/pi; 
f
%**********************繪圖****************************** 
figure('Color','w'); 
d=1:1:400; 
plot(d,M(d),d,M1(d),d,M2(d),d,M3(d),d,M4(d)); 
grid on; 
axis([0,400,0,1.1]); 
title('schmidl algorithm'); 
xlabel('Time (sample)'); 
ylabel('Timing Metric'); 
legend('no noise','SNR=1dB','SNR=5dB','SNR=10dB','SNR=10dB,f=0.45');