Split Radix FFT

DIT & DIF

通常我们见到的基-2 FFT 都是按时间抽取的 FFT，称为 Decimation-in-Time (DIT)。

下面介绍一种按频率进行抽取的基-2 FFT，称为 Decimation-in-Frequency (DIF)

X2k=∑n=0N−1exp⁡(−i2πN2kn)xn=∑n=0N/2−1exp⁡(−i2πN2kn)xn+∑n=N/2N−1exp⁡(−i2πN2kn)xn=∑n=0N/2−1exp⁡(−i2πN2kn)xn+exp⁡(−i2πN2k(n+N/2))xn+N/2=∑n=0N/2−1exp⁡(−i2πN2kn)(xn+xn+N/2)\begin{aligned} X_{2k}&=\sum\limits_{n=0}^{N-1}\exp(-i\frac{2\pi}{N}2kn)x_n\\ &=\sum\limits_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}2kn)x_n+\sum\limits_{n=N/2}^{N-1}\exp(-i\frac{2\pi}{N}2kn)x_n\\ &=\sum\limits_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}2kn)x_n+\exp\left(-i\frac{2\pi}{N}2k(n+N/2)\right)x_{n+N/2}\\ &=\sum\limits_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}2kn)(x_n+x_{n+N/2})\\ \end{aligned}

X2k+1=∑n=0N−1exp⁡(−i2πN(2k+1)n)xn=∑n=0N/2−1exp⁡(−i2πN(2k+1)n)xn+∑n=N/2N−1exp⁡(−i2πN(2k+1)n)xn=∑n=0N/2−1exp⁡(−i2πN2kn)exp⁡(−i2πNn)xn+exp⁡(−i2πN(2k+1)(n+N/2))xn+N/2=∑n=0N/2−1exp⁡(−i2πN2kn)exp⁡(−i2πNn)(xn−xn+N/2)\begin{aligned} X_{2k+1}&=\sum\limits_{n=0}^{N-1}\exp(-i\frac{2\pi}{N}\left(2k+1\right)n)x_n\\ &=\sum\limits_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}\left(2k+1\right)n)x_n+\sum\limits_{n=N/2}^{N-1}\exp(-i\frac{2\pi}{N}\left(2k+1\right)n)x_n\\ &=\sum\limits_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}2kn)\exp(-i\frac{2\pi}{N}n)x_n+\exp\left(-i\frac{2\pi}{N}\left(2k+1\right)(n+N/2)\right)x_{n+N/2}\\ &=\sum\limits_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}2kn)\exp(-i\frac{2\pi}{N}n)(x_n-x_{n+N/2}) \end{aligned}

于是可以将 {xn+xn+N/2}n=0N−1,{exp⁡(−i2πNn)(xn−xn+N/2)}n=0N−1\{x_n+x_{n+N/2}\}_{n=0}^{N-1},\{\exp(-i\frac{2\pi}{N}n)(x_n-x_{n+N/2})\}_{n=0}^{N-1} 看作两个新序列，对他们分别做 DFT。DIF 的步骤可以总结为

(a0,a1,…,aN/2−1)=(x0+xN/2,x1+x1+N/2,…,xN/2+xN−1)(b0,b1,…,bN/2−1)=(x0−xN/2,x1−x1+N/2,…,xN/2−xN−1)⊙(wN0,wN1,…,wNN/2−1)(X0,X2,…,XN−2)=FFT(a0,a1,…,aN/2−1)(X1,X3,…,XN−1)=FFT(b0,b1,…,bN/2−1)\begin{aligned} &(a_0,a_1,\dots,a_{N/2-1})=(x_0+x_{N/2},x_1+x_{1+N/2},\dots,x_{N/2}+x_{N-1})\\ &(b_0,b_1,\dots,b_{N/2-1})=(x_0-x_{N/2},x_1-x_{1+N/2},\dots,x_{N/2}-x_{N-1})\odot(w_N^0,w_N^1,\dots,w_N^{N/2-1})\\ &(X_0,X_2,\dots,X_{N-2})=\text{FFT}(a_0,a_1,\dots,a_{N/2-1})\\ &(X_1,X_3,\dots,X_{N-1})=\text{FFT}(b_0,b_1,\dots,b_{N/2-1}) \end{aligned}

DIF 与 DIT 区别在于 DIF 是先做加法和乘法，再做 FFT。而 DIT 正好相反。

Split Radix

分裂基算法是基于 DIF 推导的，实际上是基-2和基-4的混合，对偶数下标的作基-2 FFT，对奇数下标的作基-4 FFT。具体如下：

X2k=∑n=0N/2−1exp⁡(−i2πN2kn)(xn+xn+N/2)x4k+1=∑n=0N/4−1exp⁡(−i2πN4k)exp⁡(−i2πNn)(xn−xn+N/2−i(xn+N/4−xn+3N/4))x4k+3=∑n=0N/4−1exp⁡(−i2πN4k)exp⁡(−i2πN3n)(xn−xn+N/2+i(xn+N/4−xn+3N/4))\begin{aligned} X_{2k}&=\sum_{n=0}^{N/2-1}\exp(-i\frac{2\pi}{N}2kn)(x_n+x_{n+N/2})\\ x_{4k+1}&=\sum_{n=0}^{N/4-1}\exp(-i\frac{2\pi}{N}4k)\exp(-i\frac{2\pi}{N}n)(x_n-x_{n+N/2}-i(x_{n+N/4}-x_{n+3N/4}))\\ x_{4k+3}&=\sum_{n=0}^{N/4-1}\exp(-i\frac{2\pi}{N}4k)\exp(-i\frac{2\pi}{N}3n)(x_n-x_{n+N/2}+i(x_{n+N/4}-x_{n+3N/4})) \end{aligned}

这样，长度为 NN 的 DFT 可以分解为一个长度为 N/2N/2 和两个长度为 N/4N/4 的 DFT。分裂基算法相对于基-2的 FFT 计算量更小，原因是：可以观察到，偶数项前不需要乘旋转因子，这一部分计算量已经最小；而对于奇数项，由于计算机进行复数运算时是将实部和虚部分别进行计算，与 exp⁡(iπk/2),k∈Z\exp(i\pi k/2),k\in\Z 这类复数的乘积是可以简化的：exp⁡(iπk),k∈Z\exp(i\pi k),k\in\Z 不需要进行加法和乘法运算，只需要将符号或者实部虚部互换；exp⁡(iπ(k+12))\exp(i\pi(k+\frac12) ) 由于虚部和实部的绝对值相等，所以也只需进行两次加法运算和乘法运算；而其他复数则需要进行四次乘法和两次加法。

所以对于基-2 的 FFT，可以得到递推公式

Add(n)={4n=216n=42Add(n/2)+2n+2(n/2−2)=2Add(n/2)+3n−4n≥8\text{Add}(n)= \begin{cases} 4&n=2\\ 16 & n=4\\ 2\text{Add}(n/2)+2n+2(n/2-2)=2\text{Add}(n/2)+3n-4&n\ge8\\ \end{cases}\\

Mul(n)={0n=2,42Mul(n/2)+4(n/2−4)+2⋅2=2Mul(n/2)+2n−12n≥8\text{Mul}(n)= \begin{cases} 0&n=2,4\\ 2\text{Mul}(n/2)+4(n/2-4)+2\cdot2=2\text{Mul}(n/2)+2n-12&n\ge 8 \end{cases}

对于分裂基 FFT，可以得到递推公式

Add(n)={4n=216n=4Add(n/2)+2Add(n/4)+3n+2(n/2−2)n≥8\text{Add}(n)= \begin{cases} 4 & n=2\\ 16 & n=4\\ \text{Add}(n/2)+2\text{Add}(n/4)+3n+2(n/2-2)& n\ge8 \end{cases}

Mul(n)={0n=2,4Mul(n/2)+2Mul(n/4)+4(n/2−4)+4n≥8\text{Mul}(n)= \begin{cases} 0&n=2,4\\ \text{Mul}(n/2)+2\text{Mul}(n/4)+4(n/2-4)+4& n\ge8 \end{cases}

通过递推公式可以证明分裂基 FFT 的计算量，不管是加法还是减法，都更少。原因是它更好地利用了较小的 nn，同时也更好地利用了单位根的性质。

证明：

乘法是显然的，因为后面加的一项是一样的，而由计算复杂度为 Θ(nlog⁡n)\Theta(n\log n)，所以 Mul(n)≥2Mul(n/2)\text{Mul}(n)\ge2\text{Mul}(n/2)，前一项也更小

令 Add1\text{Add}_1 表示基-2 FFT，Add2\text{Add}_2 表示分裂基 FFT，运用数学归纳法，对于 n=2,4n=2,4 成立，假设对于 2,4,…,n/22,4,\dots,n/2 都成立，则有

Add1(n)−Add2(n)≥Add1(n/2)−2Add1(n/4)−n=32n−4−n=12n−4≥0\text{Add}_1(n)-\text{Add}_2(n)\ge\text{Add}_1(n/2)-2\text{Add}_1(n/4)-n=\frac32n-4-n=\frac12n-4\ge0

得证

以下展示一些 nn 的运算量

加法运算数：

n Radix-2 Split Radix
8 52 52
16 148 144
32 388 372
64 964 912

乘法运算数：

n Radix-2 Split Radix
8 4 4
16 28 24
32 108 84
64 332 248
FFT 数值算法

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！