快速傅里叶变换

定义和速度

x0, ...., xN-1复数。DFT由下式定义

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{i2\pi k{\frac {n}{N}}}}\qquad k=0,\dots ,N-1.}$

一般的简化理论

${\displaystyle \mathbf {S} ={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{m}\end{bmatrix}}{\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}}$

${\displaystyle {\begin{bmatrix}a_{1}&a_{2}&\cdots &a_{m}\end{bmatrix}}^{T}}$${\displaystyle M_{0}}$个相异的非平凡值（${\displaystyle a_{m}\neq \pm 2^{k},a_{m}\neq \pm 2^{k}a_{n}}$ where ${\displaystyle m\neq n}$）

${\displaystyle {\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}}$${\displaystyle N_{0}}$个相异的非平凡值

S共需要${\displaystyle M_{0}*N_{0}}$个乘法。

${\displaystyle {\begin{bmatrix}Z[1]\\Z[2]\\\vdots \\Z[N]\end{bmatrix}}=\mathbf {S} {\begin{bmatrix}X[1]\\X[2]\\\vdots \\X[N]\end{bmatrix}}={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{m}\end{bmatrix}}{\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}{\begin{bmatrix}X[1]\\X[2]\\\vdots \\X[N]\end{bmatrix}}}$

Step 1：${\displaystyle Z_{a}=b_{1}X[1]+b_{2}X[2]+\cdots +b_{n}X[N]}$

Step 2：${\displaystyle Z[1]=a_{1}Z_{a},Z[2]=a_{2}Z_{a},\cdots ,Z[N]=a_{m}Z_{a}}$

${\displaystyle \mathbf {S} ={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{m}\end{bmatrix}}{\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}+\mathbf {S} _{1}}$

${\displaystyle S_{1}}$也是一个M*N的矩阵。

${\displaystyle S_{1}}$${\displaystyle P_{1}}$个值不等于0，则${\displaystyle \mathbf {S} }$的乘法量上限为${\displaystyle M_{0}+N_{0}+P_{1}}$

快速傅里叶变换乘法量的计算

${\displaystyle \mathbf {P_{k}} }$点DFT的乘法量为${\displaystyle \mathbf {B_{k}} }$，则${\displaystyle \mathbf {N} }$点DFT的乘法量为：

${\displaystyle {\frac {N}{P_{1}}}B_{1}+{\frac {N}{P_{2}}}B_{2}+\cdots \cdots +{\frac {N}{P_{k}}}B_{k}}$

${\displaystyle \mathbf {N_{1}=P^{a}} }$点的DFT需要的乘法量为${\displaystyle \mathbf {B_{1}} }$

${\displaystyle n_{1}\times n_{2}}$当中（${\displaystyle n_{1}=0,1,\cdots ,N_{1}-1,\quad n_{2}=0,1,\cdots ,N_{2}-1}$

${\displaystyle D_{1}}$个值不为${\displaystyle {\frac {N}{12}}}$${\displaystyle {\frac {N}{8}}}$的倍数，

${\displaystyle D_{2}}$个值为${\displaystyle {\frac {N}{12}}}$${\displaystyle {\frac {N}{8}}}$的倍数，但不为${\displaystyle {\frac {N}{4}}}$的倍数，

${\displaystyle \mathbf {N_{2}B_{1}+N_{1}B_{2}+3D_{1}+2D_{2}} }$

库利-图基算法

设计思想

${\displaystyle W_{N}}$的性质：

1. 周期性${\displaystyle W_{N}}$具有周期${\displaystyle N}$，即${\displaystyle W_{N}^{k+N}=W_{N}^{k}}$
2. 对称性${\displaystyle W_{N}^{k+{\frac {N}{2}}}=-W_{N}^{k}}$
3. ${\displaystyle m}$${\displaystyle N}$的约数，${\displaystyle W_{N}^{mkn}=W_{\frac {N}{m}}^{kn}}$

${\displaystyle {\begin{matrix}y_{k}=\sum _{n=2t}W_{N}^{kn}x_{n}+\sum _{n=2t+1}W_{N}^{kn}x_{n}\\=\sum _{t}W_{\frac {N}{2}}^{kt}x_{2t}+W_{N}^{k}\sum _{t}W_{\frac {N}{2}}^{kt}x_{2t+1}\\=F_{even}(k)+W_{N}^{k}F_{odd}(k)&&&&&&(i\in \mathbb {Z} )\end{matrix}}}$

${\displaystyle F_{odd}(k)}$${\displaystyle F_{even}(k)}$是两个分别关于序列${\displaystyle \left\{x_{n}\right\}_{0}^{N-1}}$奇数号和偶数号序列N/2点变换。由此式只能计算出${\displaystyle y_{k}}$的前N/2个点，对于后N/2个点，注意${\displaystyle F_{odd}(k)}$${\displaystyle F_{even}(k)}$都是周期为N/2的函数，由单位根的对称性，于是有以下变换公式：

• ${\displaystyle y_{k+{\frac {N}{2}}}=F_{even}(k)-W_{N}^{k}F_{odd}(k)}$
• ${\displaystyle y_{k}=F_{even}(k)+W_{N}^{k}F_{odd}(k)}$

算法实现

• 蝶形结网络和位反转（Bit Reversal）：

001倒置以后变成 100
000 → 000
001 → 100
010 → 010
011 → 110
100 → 001
101 → 101
110 → 011
111 → 111

其他算法

Cooley-Tukey算法之外还有其他DFT的快速算法。对于长度${\displaystyle N=N_{1}N_{2}}$${\displaystyle N_{1}}$${\displaystyle N_{2}}$互质的序列，可以采用基于中国剩余定理互素因子算法将N长序列的DFT分解为两个子序列的DFT。与Cooley-Tukey算法不同的是，互素因子算法不需要旋转因子。

Bruun以及QFT算法是不断的把DFT分成许多较小的DFT运算。（Rader-Brenner以及QFT算法是为了2的指数所设计的算法，但依然可以适用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字）。尤其是Bruun算法，把FFT看成是${\displaystyle z^{N}-1}$，并把它分解成${\displaystyle z^{M-1}}$${\displaystyle z^{2M}+az^{M}+1}$的形式。

实数或对称资料专用的算法

${\displaystyle \mathbf {X} _{N-k}}$=${\displaystyle \mathbf {X} _{k}^{*}}$

function [X,Y] = Real2FFT( x,y )

% N-point FFT is used for 2 real signals both with length N

N = length(x);

z = x+y*i ;

Z = fft(z);

X = 0.5*(Z+conj(Z(1+mod(0:-1:1-N,N))));

Y = 0.5*(Z-conj(Z(1+mod(0:-1:1-N,N))))/i;

end


DFT可以应用在频谱分析以及做折积的运算，而在此处，不同应用可以用不同的算法来取代，列表如下：

• DCT
• DST
• DHT
• 正交基底的扩展（orthogonal basis expantion）包括正交多项式（orthogonal polynomials）以及CDMA.
• Haar转换
• 小波（wavelet）转换
• 时频分布（time-frequency distribution）

快速算法查表

N = 1~60

 N 乘法数 加法数 N 乘法数 N 乘法数 N 乘法数 1 0 0 11 46 24 28 39 182 2 0 4 12 8 25 148 40 100 3 2 12 13 52 26 104 42 124 4 0 16 14 32 27 114 44 160 5 10 34 15 40 28 64 45 170 6 4 36 16 20 30 80 48 92 7 16 72 18 32 32 72 52 208 8 4 52 20 40 33 160 54 228 9 16 72 21 62 35 150 56 156 10 20 88 22 80 36 64 60 160

N < 1000

 N 乘法数 N 乘法数 N 乘法数 N 乘法数 63 256 96 280 192 752 360 1540 64 204 104 468 204 976 420 2080 66 284 108 456 216 1020 480 2360 70 300 112 396 224 1016 504 2300 72 164 120 380 240 940 512 3180 80 260 128 560 252 1024 560 3100 81 480 144 436 256 1308 672 3496 84 248 160 680 288 1160 720 3620 88 412 168 580 312 1324 784 4412 90 340 180 680 336 1412 840 4580

N > 1000

 N 乘法数 N 乘法数 N 乘法数 N 乘法数 1008 5356 1440 8680 2520 16540 4032 29488 1024 7436 1680 10420 2688 19108 4096 37516 1152 7088 2016 12728 2880 20060 4368 35828 1260 7640 2048 16836 3369 24200 4608 36812 1344 8252 2304 15868 3920 29900 5040 36860

参考资料

1. ^ 杨毅明. 数字信号处理（第2版）. 北京: 机械工业出版社. 2017年: 第95页. ISBN 9787111576235.
2. ^ Charles Van Loan, Computational Frameworks for the Fast Fourier Transform (SIAM, 1992).
3. ^ Heideman, M. T.; Johnson, D. H.; Burrus, C. S. Gauss and the history of the fast Fourier transform. IEEE ASSP Magazine. 1984, 1 (4): 14–21. doi:10.1109/MASSP.1984.1162257.
4. ^ Strang, Gilbert. Wavelets. American Scientist. May–June 1994, 82 (3): 253. JSTOR 29775194.
5. ^ Dongarra, J.; Sullivan, F. Guest Editors Introduction to the top 10 algorithms. Computing in Science Engineering. January 2000, 2 (1): 22–23. ISSN 1521-9615. doi:10.1109/MCISE.2000.814652.
6. ^ Johnson and Frigo, 2007
7. ^ Frigo & Johnson, 2005