function yy=ML_cp() 
%%%Algorithm Based：J.-J. Van de Beek, M. Sandell, M. Isaksson, and P. O. B??rjesson, 
%%%“ML estimation of time and frequency offset in OFDM systems,” IEEE Trans. Signal Proc., Vol. 45, no. 7, pp. 1800-1805, July 1997. 
%%%Programmed by Xiangdong  Huang, School of  Electronics and Information Engineering,Tianjin University 
close all; clc; clear all; 
N=1024; 
L=128; 
SNR=15; 
  
disp('實際定時偏差') 
theta=87 
disp('載頻頻偏') 
e=0.25 
 
Symbols_number=6; 
Bits_number=Symbols_number*N*6; 
Bits=round(  unifrnd(0,1,1,Bits_number)     ); 
 
%%%%%%%%%%星座表 
Tab =[   -7-7*j  -7-5*j  -7-j  -7-3*j  -7+7*j  -7+5*j  -7+j  -7+3*j...   
         -5-7*j  -5-5*j  -5-j  -5-3*j  -5+7*j  -5+5*j  -5+j  -5+3*j...  
         -1-7*j  -1-5*j  -1-j  -1-3*j  -1+7*j  -1+5*j  -1+j  -1+3*j... 
         -3-7*j  -3-5*j  -3-j  -3-3*j  -3+7*j  -3+5*j  -3+j  -3+3*j... 
          7-7*j   7-5*j   7-j   7-3*j   7+7*j   7+5*j   7+j   7+3*j... 
          5-7*j   5-5*j   5-j   5-3*j   5+7*j   5+5*j   5+j   5+3*j... 
          1-7*j   1-5*j   1-j   1-3*j   1+7*j   1+5*j   1+j   1+3*j... 
          3-7*j   3-5*j   3-j   3-3*j   3+7*j   3+5*j   3+j   3+3*j ]; 
 P1=[1^2+3^2+5^2+7^2]/2/2;   
 var_s=2*P1/N; 
 SNR1=10^(SNR/10);   
 rou=SNR1/(SNR1+1); 
 lam=sqrt(var_s/(2*SNR1)); 
   
 angle_Tab=angle(Tab); 
 am_Tab=abs(Tab); 
 I_Tab=real(Tab); 
 Q_Tab=imag(Tab); 
  
 for   i=1:Symbols_number 
              input_bits=Bits((i-1)*N*6+1:i*N*6); 
              temp=ifft(QAM64_map(input_bits,Tab)); 
              s(i,:)=[temp(end-L+1:end)  temp]; 
 end 
 len=Symbols_number*(N+L); 
 s=reshape(s.',1,len); 
 k=0:len+theta-1; 
 noise=  lam*randn(1,len+theta)+j*lam*randn(1,len+theta);  
 r=[zeros(1,theta) s  ].*exp(j*2*pi*e*k/N)+noise;     %%%接收信號考慮了信道延時因素和載波頻偏因素. 
 clear noise;  clear k;  clear s; 
  
  A=zeros(1,len+theta-(2*N+L)+1); 
  ee=zeros(1,len+theta-(2*N+L)+1); 
   
  for   i=1:len+theta-(2*N+L)+1 
              r1=r(i:i+2*N+L-1); 
              for m=1:N 
                   k=m:m+L-1; 
                   gama=sum(r1(k).*conj(r1(k+N))); 
                   phi=1/2*sum( abs(r1(k)).^2+   abs(r1(k+N)).^2   );                   
              end  
               A(i)=abs(gama )-rou*phi; 
              ee(i)=-angle(gama)/(2*pi); 
  end 
   
  subplot(2,1,1); plot(0:len+theta-(2*N+L) ,A); title('似然函數曲線,峰值位置模除CP長度L即為定時位置估計');xlabel('samples') 
  subplot(2,1,2); plot(0:len+theta-(2*N+L) ,ee); hold on; title('頻偏估計曲線'); xlabel('samples') 
  plot(0:len+theta-(2*N+L) ,e*ones(1,len+theta-(2*N+L)+1),'r-.')