参考文献:
- [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.
- [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.
- [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.
- [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.
- 数字信号处理——Chirp Z变换_chirpz变换-CSDN博客
- Bluestein算法简要介绍-CSDN博客
- 如何理解离散傅里叶变换及Z变换 - 知乎 (zhihu.com)
- Bluestein’s FFT Algorithm (stanford.edu)
文章目录
- Z-transform
- Chirp Z-transform
- Bluestein's FFT
Z-transform
给定复数序列
x
n
,
n
∈
[
−
∞
,
+
∞
]
x_n, n \in [-\infty,+\infty]
xn,n∈[−∞,+∞],它的 Z-transform 定义为
X
(
z
)
=
∑
n
x
n
z
−
n
X(z) = \sum_n x_n z^{-n}
X(z)=n∑xnz−n
这是关于复变量
z
z
z 的连续函数。可以限制在有限个的非零的时域样本,
X
(
z
)
=
∑
n
=
0
N
−
1
x
n
z
−
n
X(z) = \sum_{n=0}^{N-1} x_nz^{-n}
X(z)=n=0∑N−1xnz−n
这个有限和对于任意的
z
≠
0
z\neq 0
z=0 都收敛。假如我们选取
N
N
N 个频域采样点
z
k
,
k
=
0
,
1
,
⋯
,
N
−
1
z_k,k=0,1,\cdots,N-1
zk,k=0,1,⋯,N−1,可以计算这些位置的函数值:
X
k
=
X
(
z
k
)
=
∑
n
=
0
N
−
1
x
n
z
k
−
n
X_k = X(z_k) = \sum_{n=0}^{N-1} x_nz_k^{-n}
Xk=X(zk)=n=0∑N−1xnzk−n
特别地,如果设置
z
k
=
e
2
k
π
i
/
N
z_k = e^{2k\pi i/N}
zk=e2kπi/N(复平面上单位圆的等分点),那么它恰好是 DFT
Chirp Z-transform
[RSR69] 研究了更一般的螺线上的 Z-transform,称之为啁啾 Z 变换。采样点的形式为:
z
k
=
A
W
−
k
,
k
=
0
,
1
,
⋯
,
M
−
1
z_k = AW^{-k}, k=0,1,\cdots,M-1
zk=AW−k,k=0,1,⋯,M−1
其中
M
M
M 是任意整数,
A
,
W
A,W
A,W 是任意复数,可写作如下的形式:
A
=
A
0
e
2
π
i
⋅
θ
0
W
=
W
0
e
2
π
i
⋅
ϕ
0
\begin{aligned} A &= A_0 e^{2\pi i \cdot \theta_0}\\ W &= W_0 e^{2\pi i \cdot \phi_0}\\ \end{aligned}
AW=A0e2πi⋅θ0=W0e2πi⋅ϕ0
其中,
- A 0 ∈ R + A_0 \in \mathbb R^+ A0∈R+ 表示起始采样点 z 0 z_0 z0 的模长
- θ 0 ∈ [ 0 , 1 ) \theta_0 \in [0,1) θ0∈[0,1) 确定了起始采样点 z 0 z_0 z0 的相位角,确切地说是 2 π θ 0 2\pi \theta_0 2πθ0
- ϕ 0 ∈ [ 0 , 1 ) \phi_0 \in [0,1) ϕ0∈[0,1) 确定了相邻采样点之间的角差距(angular spacing),确切地说是 2 π ϕ 0 2\pi \phi_0 2πϕ0
-
W
0
∈
R
+
W_0 \in \mathbb R^+
W0∈R+ 表示螺线的伸展率,
- 当 W 0 = 1 W_0=1 W0=1 时,采样点形成了一段圆弧
- 当 W 0 < 1 W_0<1 W0<1 时,螺旋线是外伸的(spiral out)
- 当 W 0 > 1 W_0>1 W0>1 时,螺旋线是内缩的(spiral in)
如图所示:
根据螺线上采样点的形式,可以写出:
X
k
=
X
(
A
W
−
k
)
=
∑
n
=
0
N
−
1
x
n
A
−
n
W
n
k
,
k
∈
[
M
]
X_k = X(AW^{-k}) = \sum_{n=0}^{N-1} x_nA^{-n}W^{nk},\,\, k \in [M]
Xk=X(AW−k)=n=0∑N−1xnA−nWnk,k∈[M]
啁啾 Z 变换不要求
N
,
M
N,M
N,M 相等,并且它们都是任意整数(包括素数),这里
N
N
N 是时域样本数,而
M
M
M 是频域样本数。
Bluestein’s FFT
[RSR69] 给出了计算 CZT 的快速算法。首先使用 Bluestein’s ingenious substitution,
n
k
=
n
2
+
k
2
−
(
k
−
n
)
2
2
nk = \frac{n^2 + k^2 - (k-n)^2}{2}
nk=2n2+k2−(k−n)2
于是可以得到:
X
k
=
W
k
2
2
⋅
∑
n
=
0
N
−
1
(
x
n
A
−
n
W
n
2
2
)
W
−
(
k
−
n
)
2
2
=
W
k
2
2
⋅
[
x
n
A
−
n
W
n
2
2
]
n
∈
Z
⊛
k
[
W
−
n
2
2
]
n
∈
Z
\begin{aligned} X_k &= W^{\frac{k^2}{2}} \cdot \sum_{n=0}^{N-1} \left(x_nA^{-n}W^{\frac{n^2}{2}}\right) W^{-\frac{(k-n)^2}{2}}\\ &= W^{\frac{k^2}{2}} \cdot \left[x_nA^{-n}W^{\frac{n^2}{2}}\right]_{n \in \mathbb Z} \circledast_k \left[W^{-\frac{n^2}{2}}\right]_{n \in \mathbb Z} \end{aligned}
Xk=W2k2⋅n=0∑N−1(xnA−nW2n2)W−2(k−n)2=W2k2⋅[xnA−nW2n2]n∈Z⊛k[W−2n2]n∈Z
其中符号
⊛
k
\circledast_k
⊛k 表示位置
k
∈
Z
k \in \mathbb Z
k∈Z 的卷积运算。因此 CZT 可以被这么计算:
-
首先根据序列 x n , n ∈ [ N ] x_n,n \in [N] xn,n∈[N],计算新的序列 y n = x n A − n W n 2 2 , n ∈ [ N ] y_n = x_nA^{-n}W^{\frac{n^2}{2}},\,\, n \in [N] yn=xnA−nW2n2,n∈[N]
-
其次使用 FFT 去计算(有限长的时域上)卷积运算。[Blu68] 使用了 [Sto66] 建议的 zero padding + cyclic convolution 策略,
-
设置 L = 2 m ≥ N + M − 1 L=2^m \ge N+M-1 L=2m≥N+M−1 是二的幂次
-
构造 L L L 长序列
f n = { y n , 0 ≤ n ≤ N − 1 0 , N ≤ n ≤ L − 1 f_n = \left\{\begin{aligned} y_n, && 0 \le n \le N-1\\ 0, && N \le n \le L-1 \end{aligned}\right. fn={yn,0,0≤n≤N−1N≤n≤L−1 -
构造 L L L 长序列
h n = { W − n 2 2 , 0 ≤ n ≤ M − 1 W − ( n − L ) 2 2 , M ≤ n ≤ L − 1 h_n = \left\{\begin{aligned} W^{-\frac{n^2}{2}}, && 0 \le n \le M-1\\ W^{-\frac{(n-L)^2}{2}}, && M \le n \le L-1 \end{aligned}\right. hn=⎩ ⎨ ⎧W−2n2,W−2(n−L)2,0≤n≤M−1M≤n≤L−1 -
使用 radix-2 FFT 分别计算出 F k , k ∈ [ L ] F_k, k\in[L] Fk,k∈[L] 以及 H k , k ∈ [ L ] H_k, k\in[L] Hk,k∈[L](使用 L L L-th 本原单位根),计算阿达玛乘积 G k = F k ⋅ H k G_k = F_k \cdot H_k Gk=Fk⋅Hk,再使用 Inv-FFT 计算出 g k , k ∈ [ L ] g_k, k\in[L] gk,k∈[L],并截取出前 k k k 个数值
-
-
最后计算 X k = g k ⋅ W k 2 2 , k ∈ [ M ] X_k = g_k \cdot W^{\frac{k^2}{2}},\,\, k \in [M] Xk=gk⋅W2k2,k∈[M]
总的计算开销是 O ( ( N + M ) log 2 ( N + M ) ) O\big((N+M)\log_2(N+M)\big) O((N+M)log2(N+M)),但是它花费了 3 3 3 个长度 L ≥ N + M − 1 L \ge N+M-1 L≥N+M−1 的 FFT 运算,因此隐藏的常数比 FFT 大得多。好处是时域上的采样点形成了更一般的螺线,并且数量可以是任意的,它比 DFT 更加灵活。