%********************schmidl algorithm*******************
%*******************by:sk 2008-5-22**********************
%Example:
% If
% X = rand(2,3,4);
% then
% d = size(X) returns d = [2 3 4]
% [m1,m2,m3,m4] = size(X) returns m1 = 2, m2 = 3, m3 = 4, m4 = 1
% [m,n] = size(X) returns m = 2, n = 12
% m2 = size(X,2) returns m2 = 3
close all;
clear all;
clc;
%参数定义
N=256; %FFT/IFFT 变换的点数或者子载波个数（Nu=N）
Ng=N/8; %循环前缀的长度 (保护间隔的长度)
Ns=Ng+N; %包括循环前缀的符号长度
loop=1000;
%************利用查表法生成复随机序列**********************
QAMTable=[7+7i,-7+7i,-7-7i,7-7i];
buf=QAMTable(randint(N/2,1,4)+1); %加1是为了下标可能是0不合法
%*************在奇数子载波的位置插入零*********************zj:是偶数吧？
deltaf=12.45;
SNR=0:25;
z=zeros(1,length(SNR));
SNR=0:25;
for i=1:length(SNR)
c=zeros(1,loop);
g=zeros(1,loop);
for ii=1:loop
x=zeros(N/2,1);
index = 1;
for n=1:2:N/2
x(n)=buf(index);
index=index+1;
end;
%**************利用IFFT变换生成Schmidl训练符号***************
sch = ifft(x); %[A A]的形式
sch2=[sch;(-1).*sch];
%*****************添加一个空符号以及一个后缀符号*************
src = QAMTable(randint(N,1,4)+1).';
sym = ifft(src);
sig =[zeros(N,1) sch2 sym];
%**********************添加循环前缀*************************
tx =[sig(N - Ng +1:N,:);sig];
%***********************加时延********************************
recv = reshape(tx,1,size(tx,1)*size(tx,2));
for k=1:length(recv)
recv(k)=recv(k)*exp(j*2*pi*deltaf*k/length(recv));
end
%***********************经过信道***************************
%size的1表示行，2表示列，从%前向后数，超过了为1
recv1 = awgn(recv,SNR(i),'measured');
%recv2 = awgn(recv,5,'measured');
%recv3 = awgn(recv,10,'measured');
%*****************计算符号定时*****************************
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);
for d = Ns/2+1:1:2*Ns
for k=1:2
for m=0:1:N/4-1
P(d-Ns/2) = P(d-Ns/2) + conj(recv1(d+m+(k-1)*N/2))*recv1(d+N/4+(k-1)*N/2+m);
R(d-Ns/2) = R(d-Ns/2) + power(abs(recv1(d+N/4+(k-1)*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);
end
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);
[a b]=max(M);
c(ii)=b;
end
for iii=1:loop
g(iii)=power((c(iii)-177),2);
end
z(i)=min(g)
end
figure('Color','w');
SNR=0:5:25;
semilogy(SNR,z(SNR+1),'-g*')
axis([0,25,0.00001,1000]);
xlabel('SNR(dB)');;
ylabel('Mean of Timing Offset Estimation(sample');
grid on
%**********************绘图******************************
%d=1:1:400;
%plot(d,M(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');
