快速傅里叶变换

快速傅里叶变换（英語：Fast Fourier Transform, FFT），是离散傅里叶变换的快速算法，也可用于计算离散傅里叶变换的逆变换。快速傅里叶变换有广泛的应用，如数字信号处理、计算大整数乘法、求解偏微分方程等等。本条目只描述各种快速算法。

对于复数序列$x_{0},\ x_{1},\ ...,\ x_{n-1}$，离散傅里叶变换公式为：

$y_{j}=\sum _{k=0}^{n-1}e^{-{2\pi \imath \over n}jk}x_{k}\qquad j=0,1,\dots ,n-1.$

直接变换的计算复杂度是${\mathcal {O}}(n^{2})$（参见大O符号）。快速傅里叶变换可以计算出与直接计算相同的结果，但只需要${\mathcal {O}}(n\log n)$的计算复杂度。通常，快速算法要求n能被因数分解，但不是所有的快速傅里叶变换都要求n是合数，对于所有的整数n，都存在复杂度为${\mathcal {O}}(n\log n)$的快速算法。

除了指数的符号相反、并多了一个1/n的因子，离散傅里叶变换的正变换与逆变换具有相同的形式。因此所有的离散傅里叶变换的快速算法同时适用于正逆变换。

一般的簡化理論

假設一個M*N型矩阵S可分解成列向量以及行向量相乘：

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

若${\begin{bmatrix}a_{1}&a_{2}&\cdots &a_{m}\end{bmatrix}}^{T}$有$M_{0}$個相異的非平凡值($a_{m}\neq \pm 2^{k},a_{m}\neq \pm 2^{k}a_{n}$ where $m\neq n$)

${\begin{bmatrix}b_{1}&b_{2}&\cdots &b_{n}\end{bmatrix}}$有$N_{0}$個相異的非平凡值

則S共需要$M_{0}*N_{0}$個乘法。

${\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：$Z_{a}=b_{1}X[1]+b_{2}X[2]+\cdots +b_{n}X[N]$

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

簡化理論的變型：

$\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}$

$S_{1}$也是一個M*N的矩陣。

若$S_{1}$有$P_{1}$個值不等於0，則$\mathbf {S}$的乘法量上限為$M_{0}+N_{0}+P_{1}$。

快速傅立葉變換乘法量的計算

假設$N=P_{1}\times P_{2}\times \cdots \times P_{k}$，其中$P_{1},P_{2},\cdots ,P_{k}$彼此互質

$\mathbf {P_{k}}$ 點DFT的乘法量為$\mathbf {B_{k}}$，則$\mathbf {N}$點DFT的乘法量為：

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

假設$\mathbf {N=P^{c}}$，P是一個質數。

若$\mathbf {N_{1}=P^{a}}$點的DFT需要的乘法量為$\mathbf {B_{1}}$

且$n_{1}\times n_{2}$ 當中 ($n_{1}=0,1,\cdots ,N_{1}-1,\quad n_{2}=0,1,\cdots ,N_{2}-1$)

有$D_{1}$個值不為${\frac {N}{12}}$及${\frac {N}{8}}$的倍數，

有$D_{2}$個值為${\frac {N}{12}}$及${\frac {N}{8}}$的倍數，但不為${\frac {N}{4}}$的倍數，

則N點DFT的乘法量為：

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

Cooley-Tukey算法

主条目：Cooley-Tukey快速傅里叶变换算法

Cooley-Tukey算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为$N=N_{1}N_{2}$的DFT分解为长度分别为$N_{1}$和$N_{2}$的两个较短序列的DFT，以及与${\mathcal {O}}(N)$个旋转因子的复数乘法。

这种方法以及FFT的基本思路在1965年J. W. Cooley和J. W. Tukey合作发表An algorithm for the machine calculation of complex Fourier series之后开始为人所知。但后来发现，实际上这两位作者只是重新发明了高斯在1805年就已经提出的算法(此算法在历史上数次以各种形式被再次提出)。

Cooley-Tukey算法最有名的应用，是将序列长为N 的DFT分割为两个长为N/2 的子序列的DFT，因此这一应用只适用于序列长度为2的幂的DFT计算，即基2-FFT。实际上，如同高斯和Cooley与Tukey都指出的那样，Cooley-Tukey算法也可以用于序列长度N 为任意因数分解形式的DFT，即混合基FFT，而且还可以应用于其他诸如分裂基FFT等变种。尽管Cooley-Tukey算法的基本思路是采用递归的方法进行计算，大多数传统的算法实现都将显式的递归算法改写为非递归的形式。另外，因为Cooley-Tukey算法是将DFT分解为较小长度的多个DFT，因此它可以同任一种其他的DFT算法联合使用。

设计思想

下面，我们用N次单位根$W_{N}$来表示$e^{-j{\frac {2\pi }{N}}}$。

$W_{N}$的性质：

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

为了简单起见，我们下面设待变换序列长度$n=2^{r}$。 根据上面单位根的对称性，求级数$y_{k}=\sum _{n=0}^{N-1}W_{N}^{kn}x_{n}$时，可以将求和区间分为两部分：

${\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}}$

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

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

这样，一个N点变换就分解成了两个N/2点变换。照这样可继续分解下去。这就是库利-图基快速傅里叶变换算法的基本原理。根据主定理不难分析出此时算法的时间复杂度为$O(N\log N)$

其他算法

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

Rader-Brenner 算法是类似于 Cooley-Tukey 算法，但是采用的旋转因子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radix variant of Cooley-Tukey所取代，与Rader-Brenner演算法相比，有一样多的乘法量，却有较少的加法量，且不牺牲数值的准确性。

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

另一个从多项式观点的快速傅立叶变换法是Winograd 算法。此演算法把$z^{N}-1$分解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。这样只需要很少的乘法量(如果有需要的话)，所以winograd是可以得到最少乘法量的快速傅立叶演算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd演算法让DFT可以用$2^{k}$点的DFT来简化，但减少乘法量的同时，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。

Rader演算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋摺积来表示原本的DFT，如此就可利用摺积用一对基本的FFT来计算DFT。另一个prime-size的FFT演算法为chirp-Z演算法。此法也是将DFT用摺积来表示，此法与Rader演算法相比，能运用在更一般的转换上，其转换的基础为Z转换(Rabiner et al., 1969)。

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

在许多的运用当中，要进行DFT的资料是纯实数，如此一来经过DFT的结果会满足对称性：

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

而有一些演算法是专门为这种情形设计的(e.g. Sorensen, 1987)。另一些则是利用旧有的演算法(e.g. Cooley-Tukey)，删去一些不必要的演算步骤，如此省下了记忆体的使用，也省下了时间。另一方面，也可以把一个偶数长度且纯实数的DFT，用长度为原本一半的复数型态DFT来表示(实数项为原本纯实数资料的偶数项，虚数项则为奇数项)。

一度人们认为，用离散哈特利转换(Discrete Hartley Transform)来处理纯实数的DFT会更有效率，但接着人们发现，对于同样点数的纯实数DFT，经过设计的FFT，可以比DHT省下更多的运算。Bruun演算法是第一个试着从减少实数输入的DFT运算量的演算法，但此法并没有成为人们普遍使用的方法。

对于实数输入，且输入为偶对称或奇对称的情形，可以更进一步的省下时间以及记忆体，此时DFT可以用离散余弦转换或离散正弦转换来代替(Discrete cosine/sine transforms)。由于DCT/DST也可以设计出FFT的演算法，故在此种情形下，此方法取代了对DFT设计的FFT演算法。

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

用来做频谱分析的情况下，DFT可用下列的演算法代替：

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

用来做摺积的情况下，DFT可用下列的演算法代替：

• DCT.
• DST.
• DHT.
• 直接做摺积(direct computing)
• 分段式DFT摺积(sectioned DFT convolution)
• 威諾格拉德快速傅立葉變換演算法
• Walsh(Hadamard)转换
• 数论转换

複雜度以及運算量的極限

長久以來，人們對於求出快速傅立葉變換的複雜度下限以及需要多少的運算量感到很有興趣，而實際上也還有許多問題需要解決。即使是用較簡單的情形，即$2^{k}$點的DFT，也還沒能夠嚴謹的證明出FFT至少需要$\Omega (NlogN)$(比NlogN大)的運算量，目前也沒有發現複雜度更低的演算法。通常數學運算量的多寡會是運算效率好壞最主要的因素，但在現實中，有許多因素也會有很大的影響，如快取記憶體以及CPU均有很大的影響。

在1978年，Winograd率先導出一個較嚴謹的FFT所需乘法量的下限：$\Theta (N)$。當$N=2^{k}$時，DFT只需要$4N-2\log _{2}^{2}N-2\log _{2}N-4$ 次無理實數的乘法即可以計算出來。更詳盡，且也能趨近此下限的演算法也一一被提出(Heideman & Burrus, 1986; Duhamel, 1990)。很可惜的是，這些演算法，都需要很大量的加法計算，目前的硬體無法克服這個問題。

對於所需加法量的數目，雖然我們可以在某些受限制的假設下，推得其下限，但目前並沒有一個精確的下限被推導出來。1973年，Morgenstern在乘法常數趨近巨大的情形下(對大部分的FFT演算法為真，但不是全部)推導出加法量的下限：$\Omega \left(N\log N\right)$。Pan(1986)在假設FFT演算法的不同步的情形有其極限下證明出加法量的下限$\Omega (NlogN)$，但一般來說，此假設相當的不明確。長度為$N=2^{k}$的情形下，在某些假設下，Papadimitriou(1979)提出使用Cooley-Tukey演算法所需的複數加法量$N\log _{2}N$是最少的。到目前為止，在長度為$N=2^{k}$情況，還沒有任何FFT的演算法可以讓複數的加法量比$N\log _{2}N$還少。

還有一個問題是如何把乘法量與加法量的總和最小化，有時候稱作"演算複雜度"(在這裡考慮的是實際的運算量，而不是漸近複雜度)。同樣的，沒有一個嚴謹下限被證明出來。從1968年開始，$N=2^{k}$點DFT而言，split-radix FFT演算法需要最少的運算量，在$N>1$的情形下，其需要$4N\log _{2}N-6N+8$ 個乘法運算以及加法運算。最近有人導出更低的運算量：${\frac {34}{9}}N\log _{2}N$。(Johnson and Frigo, 2007; Lundy and Van Buskirk, 2007)

大多數嘗試要降低或者證明FFT複雜度下限的人都把焦點放在複數資料輸入的情況，因其為最簡單的情形。但是，複數資料輸入的FFT演算法，與實數資料輸入的FFT演算法，離散餘旋轉換(DCT)，離散哈特列轉換(DHT)，以及其他的演算法，均有很大的關連性。故任何一個演算法，在複雜度上有任何的改善的話，其他的演算法複雜度也會馬上獲得改善(Duhamel & Vetterli, 1990)。

参考资料

• N. Brenner and C. Rader, 1976, A New Principle for Fast Fourier Transformation, IEEE Acoustics, Speech & Signal Processing 24: 264-266.
• Cooley, James W., and John W. Tukey, 1965, "An algorithm for the machine calculation of complex Fourier series," Math. Comput. 19: 297–301.
• Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein, 2001. Introduction to Algorithms, 2nd. ed. MIT Press and McGraw-Hill. ISBN 0-262-03293-7. Especially chapter 30, "Polynomials and the FFT."
• Pierre Duhamel, 1990, , IEEE Trans. Acoust. Speech. Sig. Proc. 38: 1504-151.
• ------- and M. Vetterli, 1990, , Signal Processing 19: 259–299.
• A. Edelman, P. McCorquodale, and S. Toledo, 1999, , SIAM J. Sci. Computing 20: 1094–1114.
• Funda Ergün, 1995, , Proc. 27th ACM Symposium on the Theory of Computing: 407–416.
• M. Frigo and S. G. Johnson, 2005, "The Design and Implementation of FFTW3," Proceedings of the IEEE 93: 216–231.
• Carl Friedrich Gauss, 1866. "Nachlass: Theoria interpolationis methodo nova tractata," Werke band 3, 265–327. Göttingen: Königliche Gesellschaft der Wissenschaften.
• W. M. Gentleman and G. Sande, 1966, "Fast Fourier transforms—for fun and profit," Proc. AFIPS 29: 563–578.
• H. Guo and C. S. Burrus, 1996, , Proc. SPIE Intl. Soc. Opt. Eng. 2825: 250–259.
• ------- and G. A. Sitton, 1994, , Proc. IEEE Conf. Acoust. Speech and Sig. Processing (ICASSP) 3: 445–448.
• Michael T. Heideman and C. Sidney Burrus, 1986, On the number of multiplications necessary to compute a length-$2^{n}$ DFT, IEEE Trans. Acoust. Speech. Sig. Proc. 34: 91-95.
• -------- and D. H. Johnson, 1984, Gauss and the history of the fast Fourier transform, IEEE ASSP Magazine 1: 14–21.
• S. G. Johnson and M. Frigo, 2007. "A modified split-radix FFT with fewer arithmetic operations," IEEE Trans. Signal Processing 55 (1): 111–119.
• T. Lundy and J. Van Buskirk, 2007. "A new matrix approach to real FFTs and convolutions of length 2^k," Computing 80 (1): 23-45.
• Jacques Morgenstern, 1973, , J. ACM 20: 305-306.
• M. J. Mohlenkamp, 1999, "A fast transform for spherical harmonics", J. Fourier Anal. Appl. 5, 159–184. (preprint)
• H. J. Nussbaumer, 1977, , Electronics Lett. 13: 386-387.
• V. Pan, 1986, , Information Proc. Lett. 22: 11-14.
• Christos H. Papadimitriou, 1979, , J. ACM 26: 95-102.
• D. Potts, G. Steidl, and M. Tasche, 2001. "Fast Fourier transforms for nonequispaced data: A tutorial", in: J.J. Benedetto and P. Ferreira (Eds.), Modern Sampling Theory: Mathematics and Applications (Birkhauser).
• Vladimir Rokhlin and Mark Tygert, 2006, "Fast algorithms for spherical harmonic expansions," SIAM J. Sci. Computing 27 (6): 1903-1928.
• James C. Schatzman, 1996, Accuracy of the discrete Fourier transform and the fast Fourier transform, SIAM J. Sci. Comput. 17: 1150–1166.
• O. V. Shentov, S. K. Mitra, U. Heute, and A. N. Hossen, 1995, , Signal Processing 41: 261–277.
• H. V. Sorensen, D. L. Jones, M. T. Heideman, and C. S. Burrus, 1987, Real-valued fast Fourier transform algorithms, IEEE Trans. Acoust. Speech Sig. Processing ASSP-35: 849–863. See also Corrections to "Real-valued fast Fourier transform algorithms"
• Peter D. Welch, 1969, A fixed-point fast Fourier transform error analysis, IEEE Trans. Audio Electroacoustics 17: 151–157.
• S. Winograd, 1978, On computing the discrete Fourier transform, Math. Computation 32: 175-199.
• Jian-Jiun Ding, Advanced Digital Signal Processing class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2007.
• E. O. Brigham, The　Fast Fourier Transform,Prentice Hall,Englewood Cliffs,New Jersey,1974.
• E.O.布赖姆著,柳群译,快速富里叶变换,上海科学技术出版社,1979.

参阅

• 离散傅里叶变换
• 并行快速傅里叶变换