Chirp-Z变换（Chirp-Z Transform，CZT）

采用FFT算法可以很快地计算出全部DFT值，即Z变换在单位圆上的全部等间隔采样值。
在实际情况中，并不需要对整个单位圆的频谱进行分析，例如，对于窄带信号，往往只需要对信号所在的一段频带进行分析，即可在所关心的这段频带内进行密集的采样，而对这个频带以外的部分可以完全不管。
Z变换的螺旋采样，它沿Z平面上的一段螺线进行等分角的采样，这些采样点可以表示为
$$z_k=A{W}^{-k},k=0,1,\cdots ,M-1$$
其中，$M$为采样点的总数，$A$为起始点位置，可以用半径$A_0$及相角${\theta }_0$表示为$A=A_0{e}^{j{\theta }_0}$；
参数$W$表示为$W=W_0{e}^{-j{\varphi }_0}$，$W_0$为螺线的伸展率，$W_0>1$，螺线内缩，$W_0<1$，螺线外伸;
${\varphi }_0$为螺线上采样点之间的等分角。

当$M=N$、$A=1$、$W=e^{-j\frac{2\pi }{N}}$时，$z_k$就等间隔地分布在单位圆上，这时CZT退化DFT。
假设$x(n)$是长度为$N$的有限长序列，则其Z变换在采样点$z_k$上的值$$X(z_k)=\sum _{n=0}^{N-1}x(n)z_k^{-n},k=0,1,1,\cdots ,M-1$$
为减少计算量，将上式运算转换为卷积形式，从而采用FFT进行计算。

算法原理

将$z_k=A{W}^{-k}$代入$X(z_k)={\sum }_{n=0}^{N-1}x(n)z_k^{-n}$可得
$$X(z_k)=\sum _{n=0}^{N-1}x(n)A^{-n}{W}^{nk}$$
将$nk$替换为$\frac{1}{2}[k^2+n^2-(k-n{)}^2]$，则
$$X(z_k)=\sum _{n=0}^{N-1}x(n)A^{-n}{W}^{\frac{1}{2}[k^2+n^2-(k-n{)}^2]}=W^{\frac{k^2}{2}}\sum _{n=0}^{N-1}x(n)A^{-n}{W}^{\frac{n^2}{2}}{W}^{-\frac{(k-n{)}^2}{2}}$$
定义$g(n)=x(n)A^{-n}{W}^{\frac{n^2}{2}},n=0,1,2,\cdots ,N-1$和$h(n)=W^{\frac{-n^2}{2}}$，则有
$$g(k)\ast h(k)=\sum _{n=0}^{N-1}g(n)h(k-n)=\sum _{n=0}^{N-1}x(n)A^{-n}{W}^{\frac{n^2}{2}}{W}^{-\frac{(k-n{)}^2}{2}},k=0,1,\cdots ,M-1$$
则有
$$X(z_k)=[g(k)\ast h(k)]W^{k^2}2,k=0,1,\cdots ,M-1$$
算法流程图如下：

实现步骤

1. 选择一个最小整数$L$，使其满足$L\ge N+M-1$，同时$L=2^m$；
2. 求$h(n)$的主值序列$\hat{h}(n)$，并计算DFT；
   $$\hat{h}(n)=\{\begin{array}{ll}{W}^{-\frac{n^2}{2}}& 0\le n\le M-1\\ \text{任意值}& N\le n\le L-1\\ W^{-\frac{(n-L{)}^2}{2}}& L-N+1\le n\le L-1\end{array}$$
   $$H(k)=DFT[\hat{h}(n)],L点$$
3. 对$x(n)$加权、补零，并计算DFT；
   $$g(n)=\{\begin{array}{ll}x(n)A^{-n}{W}^{\frac{n^2}{2}}& 0\le n\le N-1\\ 0& N\le nL-1\end{array}$$
   $$G(k)=DFT[g(n)],L点$$
4. $Y(k)=G(k)H(k),L点$；
5. $y(n)=IDFT[Y(k)],L点$；
6. $X(z_k)=W^{\frac{k^2}{2}}y(k),0\le k\le M-1$。

上述步骤实现程序可见Matlab的czt函数内部程序。

仿真分析

此处使用函数czt实现Chirp-Z变换，并将结果与DFT和采样序列插0后序列的DFT进行对比。

clc;clear;close all;
N = 8192;
f1 = 100;
f2 = 101;
fs = 8000;
Ts = 1/fs; 
ts = (1:N)*Ts;
x = cos(2*pi*f1*ts) + cos(2*pi*f2*ts) + 0.5*randn(1,N); 
y_DFT = abs(fft(x)); %%DFT
w = exp(-1i*2*pi*(150-50)/(N*fs));
a = exp(1i*2*pi*50/fs);
y_CZT = abs(czt(x,N,w,a));%%CZT

fn = (0:N-1)/N;
fy = fs*fn;
fz = (150-50)*fn + 50;
fyy = fs*(0:2*N-1)/(2*N);
xx = [x zeros(1,N)];
yy_DFT = abs(fft(xx)); %%插0 DFT
plot(yy_DFT);
plot(fy,20*log10(y_DFT), fz,20*log10(y_CZT), fyy,20*log10(yy_DFT));
xlim([80 120]);
legend('DFT','CZT','插0 DFT');