mczt.m

对低频段的带通信号进行高精度的时延估计

%期刊中的MCZT算法
function Xk = MCZT(xn,N,N1)
%MCZT函数表示N1等分2*pi时，对N点时间序列xn计算点N点的chirp-z变换。
%Xk表示N点的变换输出序列
%IMCZT是所要变换序列共轭的MCZT的共轭，即IMCZT（Xk）=conj[MCZT(conj(Xk))]
%MCZT的计算过程如下:
%  1-----计算权函数h0(n)和W(n) n=0-（N-1)
n=0:N-1;
h0=exp(j*pi*n.*n/N1);
Wn=conj(h0);
%Hw=hamming(N);     xn=xn.*Hw';  %加窗反而会使性能变差
f1=xn.*Wn;
%  2-----选择合适的L值，L大于或者等于2*N-1，而且又要是2的幂
%  NEXTPOW2(N) returns the first P such that 2^P >= abs(N).
Ln=nextpow2(2*N-1);
L=2^Ln;
%  3-----补零将f1扩充为L点的序列f(n)
f=[f1,zeros(1,L-N)];
%  4-----补零将h0(n)扩充为L点的序列h(n)
n2=L-N:-1:1;%计算L-N点的值
h1=exp(j*pi*n2.*n2/N1);
h=[h0,h1];
%  5-----利用FFT来计算f(n)和h(n)的圆周卷积,都是L点的运算
Fk=fft(f,L);
Hk=fft(h,L);
Gk=Fk.*Hk;
Gn=ifft(Gk,L);
%  6-----取前面的N点序列gk
gk(1:N)=Gn(1:N);
%gk(1:3105)=0;
Xk=gk.*Wn;
%Xk(1:3105)=0;