czt.m

基于FFT的极大似然载频估计算法

fc=3000;
% fk=zeros(1,length(fc));
% for i=1:length(fc)    
    fi=16e3;
    fr=fc+fi;
fs=(1.28e6)*4;
t=1/(100*fs):1/fs:513/fs;
re=cos(2*pi*fr*t);
ev=(15e3:21e3);
fv1=min(ev);
fv2=max(ev);
tell=50;
 A0=1;
 W0=1;
 sita=2*pi*fv1/fs;
 N=513;
 M=(fv2-fv1)/tell+1;
 L=1024;
 fai=2*pi*(fv2-fv1)/(fs*(M-1));
 A=A0*exp(j*sita);
 W=W0*exp(-j*fai);
% h=zeros(1,1024);
%  for n=1:L
%     if n<=M
%     h(n)=W^(-(n-1)^2/2);
%     else
%     h(n)=W^(-(n-L-1)^2/2);
%     end
% end
% 
% for s=1:N
%     y1(s)=re(s)*(A^(-s+1))*(W^((s-1)^2/2));
% end
% y=[y1,zeros(1,L-N)];
% H=fft(h,L);
% Y=fft(y,L);
% V=H.*Y;
% v=ifft(V,L);
%  k=0:L-1;
% out=W.^(k.^2/2).*v;
% out=out(1:M);
% %out=fft(re,65536);
% p=out.*conj(out);
% [c,m]=max(p);
% % fe = (m-1)/(65536*fs)- fi
% % fk=(sita+fai*(m))/(fs*2*pi)-fi
% % fe = (m-1)/(65536*fs)- fi
% % fk(r)=(sita+fai*(m-1))/(fs*2*pi)-fi;
% % fe = (m-1)/(65536*fs)- fi
n = (0:M-1);
hl = W .^ (-n .^ 2 ./ 2);
n = (L-N+1:L-1);
hltmp = W .^ (-(n-L) .^ 2 ./2);
hl = [hl,zeros(1,L-N-M+1),hltmp];

n = (0:N-1);
y = re .* A .^(-n) .* W .^(n.^2/2);
y = [y,zeros(1,L-N)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算X(k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = fft(hl);
Y = fft(y);
V = H .* Y;
v = ifft(V);
n = (0:L-1);
Xtmp = W .^ (n.^2/2) .* v;
X = Xtmp(1:M);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%输出估计值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pyy = X.* conj(X) / L;
[c,m] = max(Pyy);
% fe = (sita + (m-1)*fai)*fs / (2*pi)

fk=(sita+fai*(m-1))*fs/(2*pi)-fi
% end
% plot(fc,fk,'b:p')
% title('the performance of czt')
% xlabel('frequency offset (Hz)');
% ylabel('czt output (Hz)');