CZT & Blusetein‘s FFT

参考文献：

1. [Sto66] Stockham Jr T G. High-speed convolution and correlation[C]//Proceedings of the April 26-28, 1966, Spring joint computer conference. 1966: 229-233.
2. [Blu68] Bluestein L. A linear filtering approach to the computation of discrete Fourier transform[J]. IEEE Transactions on Audio and Electroacoustics, 1970, 18(4): 451-455.
3. [RSR69] Rabiner L R, Schafer R W, Rader C M. The chirp z‐transform algorithm and its application[J]. Bell System Technical Journal, 1969, 48(5): 1249-1292.
4. [PKG+16] Pariyal P S, Koyani D M, Gandhi D M, et al. Comparison based analysis of different FFT architectures[J]. International journal of image, graphics and signal processing, 2016, 8(6): 41.
5. 数字信号处理——Chirp Z变换_chirpz变换-CSDN博客
6. Bluestein算法简要介绍-CSDN博客
7. 如何理解离散傅里叶变换及Z变换 - 知乎 (zhihu.com)
8. Bluestein’s FFT Algorithm (stanford.edu)

Z-transform

给定复数序列 $x_n,n\in [-\infty ,+\infty ]$，它的 Z-transform 定义为
$$X(z)=\sum _n{x}_n{z}^{-n}$$
这是关于复变量 $z$ 的连续函数。可以限制在有限个的非零的时域样本，
$$X(z)=\sum _{n=0}^{N-1}{x}_n{z}^{-n}$$
这个有限和对于任意的 $z\ne 0$ 都收敛。假如我们选取 $N$ 个频域采样点 $z_k,k=0,1,\cdots ,N-1$，可以计算这些位置的函数值：
$$X_k=X(z_k)=\sum _{n=0}^{N-1}{x}_n{z}_k^{-n}$$
特别地，如果设置 $z_k=e^{2k\pi i/N}$（复平面上单位圆的等分点），那么它恰好是 DFT

Chirp Z-transform

[RSR69] 研究了更一般的螺线上的 Z-transform，称之为啁啾 Z 变换。采样点的形式为：
$$z_k=A{W}^{-k},k=0,1,\cdots ,M-1$$
其中 $M$ 是任意整数，$A,W$ 是任意复数，可写作如下的形式：
$$\begin{array}{rl}A& =A_0{e}^{2\pi i\cdot {\theta }_0}\\ W& =W_0{e}^{2\pi i\cdot {\varphi }_0}\end{array}$$
其中，

• $A_0\in {\mathbb{R}}^{+}$ 表示起始采样点 $z_0$ 的模长
• ${\theta }_0\in [0,1)$ 确定了起始采样点 $z_0$ 的相位角，确切地说是 $2\pi {\theta }_0$
• ${\varphi }_0\in [0,1)$ 确定了相邻采样点之间的角差距（angular spacing），确切地说是 $2\pi {\varphi }_0$
• $W_0\in {\mathbb{R}}^{+}$ 表示螺线的伸展率，
  • 当 $W_0=1$ 时，采样点形成了一段圆弧
  • 当 $W_0<1$ 时，螺旋线是外伸的（spiral out）
  • 当 $W_0>1$ 时，螺旋线是内缩的（spiral in）

如图所示：

根据螺线上采样点的形式，可以写出：
$$X_k=X(A{W}^{-k})=\sum _{n=0}^{N-1}{x}_n{A}^{-n}{W}^{nk},k\in [M]$$
啁啾 Z 变换不要求 $N,M$ 相等，并且它们都是任意整数（包括素数），这里 $N$ 是时域样本数，而 $M$ 是频域样本数。

Bluestein’s FFT

[RSR69] 给出了计算 CZT 的快速算法。首先使用 Bluestein’s ingenious substitution，
$$nk=\frac{n^2+k^2-(k-n{)}^2}{2}$$
于是可以得到：
$$\begin{array}{rl}{X}_k& =W^{\frac{k^2}{2}}\cdot \sum _{n=0}^{N-1}\left(x_n{A}^{-n}{W}^{\frac{n^2}{2}}\right)W^{-\frac{(k-n{)}^2}{2}}\\ & =W^{\frac{k^2}{2}}\cdot {\left[x_n{A}^{-n}{W}^{\frac{n^2}{2}}\right]}_{n\in \mathbb{Z}}{⊛}_k{\left[W^{-\frac{n^2}{2}}\right]}_{n\in \mathbb{Z}}\end{array}$$
其中符号 ${⊛}_k$ 表示位置 $k\in \mathbb{Z}$ 的卷积运算。因此 CZT 可以被这么计算：

1. 首先根据序列 $x_n,n\in [N]$，计算新的序列 $y_n=x_n{A}^{-n}{W}^{\frac{n^2}{2}},n\in [N]$

2. 其次使用 FFT 去计算（有限长的时域上）卷积运算。[Blu68] 使用了 [Sto66] 建议的 zero padding + cyclic convolution 策略，

   1. 设置 $L=2^m\ge N+M-1$ 是二的幂次

   2. 构造 $L$ 长序列
      $$f_n=\{\begin{array}{rlr}{y}_n,& & 0\le n\le N-1\\ 0,& & N\le n\le L-1\end{array}$$

   3. 构造 $L$ 长序列
      $$h_n=\{\begin{array}{rlr}{W}^{-\frac{n^2}{2}},& & 0\le n\le M-1\\ W^{-\frac{(n-L{)}^2}{2}},& & M\le n\le L-1\end{array}$$

   4. 使用 radix-2 FFT 分别计算出 $F_k,k\in [L]$ 以及 $H_k,k\in [L]$（使用 $L$-th 本原单位根），计算阿达玛乘积 $G_k=F_k\cdot H_k$，再使用 Inv-FFT 计算出 $g_k,k\in [L]$，并截取出前 $k$ 个数值

3. 最后计算 $X_k=g_k\cdot W^{\frac{k^2}{2}},k\in [M]$

总的计算开销是 $O((N+M){\log }_2(N+M))$，但是它花费了 $3$ 个长度 $L\ge N+M-1$ 的 FFT 运算，因此隐藏的常数比 FFT 大得多。好处是时域上的采样点形成了更一般的螺线，并且数量可以是任意的，它比 DFT 更加灵活。