Brstjjy

  「FFT」重定向至此。關於与此名称相似的其他条目，詳見「FFT (消歧义)」。

一个五项余弦级数的时域信号。频率为 102{displaystyle 10{sqrt {2}}} 的倍数。振幅为 1, 2, 3, 4, 5。时间步长为 0.001 s

用FFT（~40,000次运算）五项余弦级数及用显式积分（~1亿次运算）得出的DFT。时间窗口是10秒。FFT是用FFTW3（ http://www.fftw.org/ ）计算的。显式积分DFT是用 https://sourceforge.net/projects/amoreaccuratefouriertransform/ 计算的。

缩放后的五项余弦级数的DFT。注意到显式积分更细的步长大小比FFT更精确地再现了峰值（4）和频率（56.569 Hz），代价是速度慢了上千倍。

快速傅里叶变换（英语：Fast Fourier Transform, FFT），是快速计算序列的离散傅里叶变换（DFT）或其逆变换的方法^[1]。傅里叶分析将信号从原始域（通常是时间或空间）转换到頻域的表示或者逆过来转换。FFT会通过把DFT矩阵分解为稀疏（大多为零）因子之积来快速计算此类变换。^[2] 因此，它能够将计算DFT的复杂度从只用DFT定义计算需要的 O(n2){displaystyle O(n^{2})}，降低到 O(nlog⁡n){displaystyle O(nlog n)}，其中 n{displaystyle n} 为数据大小。

快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及，但早在1805年就已推导出来。^[3] 1994年美國數學家吉尔伯特·斯特朗（英语：Gilbert Strang）把FFT描述为“我们一生中最重要的数值算法”^[4]，它还被IEEE科学与工程计算期刊列入20世纪十大算法。^[5]

定义和速度

用FFT计算DFT会得到与直接用DFT定义计算相同的结果；最重要的区别是FFT更快。（由于捨入誤差的存在，许多FFT算法还会比直接运用定义求值精确很多，后面会讨论到这一点。）

令 x₀, ...., x_N-1 為复数。DFT由下式定义

Xk=∑n=0N−1xne−i2πknNk=0,…,N−1.{displaystyle X_{k}=sum _{n=0}^{N-1}x_{n}e^{-{i2pi k{frac {n}{N}}}}qquad k=0,dots ,N-1.}

直接按这个定义求值需要 O(N²) 次运算：X_k 共有 N 个输出，每个输出需要 N 项求和。直接使用DFT運算需使用N個複數乘法(4N 個實數乘法)與N-1個複數加法(4N-4個實數加法)，因此，計算使用DFT所有N點的值需要N²複數乘法與N²-N 個複數加法。FFT则是能够在 O(N log N) 次操作计算出相同结果的任何方法。更准确的说，所有已知的FFT算法都需要 O(N log N) 次运算（技术上O只标记上界），虽然还没有已知的证据证明更低的复杂度是不可能的。（Johnson and Frigo, 2007）

要说明FFT节省时间的方式，就得考虑复数相乘和相加的次数。直接计算DFT的值涉及到 N² 次复数相乘和 N(N−1) 次复数相加（可以通过削去琐碎运算（如乘以1）来节省 O(N) 次运算）。众所周知的基2库利-图基算法，N 为2的幂，可以只用 (N/2)log₂(N) 次复数乘法（再次忽略乘以1的简化）和 Nlog₂(N) 次加法就可以得到相同结果。在实际中，现代计算机通常的实际性能通常不受算术运算的速度和对复杂主体的分析主导（参见Frigo & Johnson, 2005），但是从 O(N²) 到 O(N log N) 的总体改进仍然能够体现出来。

一般的簡化理論

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

S=[a1a2⋮am][b1b2⋯bn]{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}}}

若[a1a2⋯am]T{displaystyle {begin{bmatrix}a_{1}&a_{2}&cdots &a_{m}end{bmatrix}}^{T}}有M0{displaystyle M_{0}}個相異的非平凡值（am≠±2k,am≠±2kan{displaystyle a_{m}neq pm 2^{k},a_{m}neq pm 2^{k}a_{n}} where m≠n{displaystyle mneq n}）　

[b1b2⋯bn]{displaystyle {begin{bmatrix}b_{1}&b_{2}&cdots &b_{n}end{bmatrix}}}有N0{displaystyle N_{0}}個相異的非平凡值

則S共需要M0∗N0{displaystyle M_{0}*N_{0}}個乘法。

[Z[1]Z[2]⋮Z[N]]=S[X[1]X[2]⋮X[N]]=[a1a2⋮am][b1b2⋯bn][X[1]X[2]⋮X[N]]{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：Za=b1X[1]+b2X[2]+⋯+bnX[N]{displaystyle Z_{a}=b_{1}X[1]+b_{2}X[2]+cdots +b_{n}X[N]}

Step 2：Z[1]=a1Za,Z[2]=a2Za,⋯,Z[N]=amZa{displaystyle Z[1]=a_{1}Z_{a},Z[2]=a_{2}Z_{a},cdots ,Z[N]=a_{m}Z_{a}}

簡化理論的變型：

S=[a1a2⋮am][b1b2⋯bn]+S1{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}}

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

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

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

假設N=P1×P2×⋯×Pk{displaystyle N=P_{1}times P_{2}times cdots times P_{k}}，其中P1,P2,⋯,Pk{displaystyle P_{1},P_{2},cdots ,P_{k}}彼此互質

Pk{displaystyle mathbf {P_{k}} }點DFT的乘法量為Bk{displaystyle mathbf {B_{k}} }，則N{displaystyle mathbf {N} }點DFT的乘法量為：

NP1B1+NP2B2+⋯⋯+NPkBk{displaystyle {frac {N}{P_{1}}}B_{1}+{frac {N}{P_{2}}}B_{2}+cdots cdots +{frac {N}{P_{k}}}B_{k}}

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

若N1=Pa{displaystyle mathbf {N_{1}=P^{a}} }點的DFT需要的乘法量為B1{displaystyle mathbf {B_{1}} }

且n1×n2{displaystyle n_{1}times n_{2}}當中（n1=0,1,⋯,N1−1,n2=0,1,⋯,N2−1{displaystyle n_{1}=0,1,cdots ,N_{1}-1,quad n_{2}=0,1,cdots ,N_{2}-1}）

有D1{displaystyle D_{1}}個值不為N12{displaystyle {frac {N}{12}}}及N8{displaystyle {frac {N}{8}}}的倍數，

有D2{displaystyle D_{2}}個值為N12{displaystyle {frac {N}{12}}}及N8{displaystyle {frac {N}{8}}}的倍數，但不為N4{displaystyle {frac {N}{4}}}的倍數，

則N點DFT的乘法量為：

N2B1+N1B2+3D1+2D2{displaystyle mathbf {N_{2}B_{1}+N_{1}B_{2}+3D_{1}+2D_{2}} }

库利-图基算法

库利-图基算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为N=N1N2{displaystyle N=N_{1}N_{2}}的离散傅里叶变换分解为长度为N1{displaystyle N_{1}}的N2{displaystyle N_{2}}个较短序列的离散傅里叶变换，以及与O(N){displaystyle mathrm {O} (N)}个旋转因子的复数乘法。

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

库利-图基算法最有名的应用，是将序列长为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次单位根WN{displaystyle W_{N}}来表示e−j2πN{displaystyle e^{-j{frac {2pi }{N}}}}。

WN{displaystyle W_{N}}的性质：

1. 
   周期性，WN{displaystyle W_{N}}具有周期N{displaystyle N}，即WNk+N=WNk{displaystyle W_{N}^{k+N}=W_{N}^{k}}
   


2. 
   对称性：WNk+N2=−WNk{displaystyle W_{N}^{k+{frac {N}{2}}}=-W_{N}^{k}}。


3. 若m{displaystyle m}是N{displaystyle N}的约数，WNmkn=WNmkn{displaystyle W_{N}^{mkn}=W_{frac {N}{m}}^{kn}}
   

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

yk=∑n=2tWNknxn+∑n=2t+1WNknxn=∑tWN2ktx2t+WNk∑tWN2ktx2t+1=Feven(k)+WNkFodd(k)(i∈Z){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)&&&&&&(iin mathbb {Z} )end{matrix}}}

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

• yk+N2=Feven(k)−WNkFodd(k){displaystyle y_{k+{frac {N}{2}}}=F_{even}(k)-W_{N}^{k}F_{odd}(k)}


• 
  yk=Feven(k)+WNkFodd(k){displaystyle y_{k}=F_{even}(k)+W_{N}^{k}F_{odd}(k)}。

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

算法实现

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

首先将n=2N{displaystyle n=2^{N}}个输入点列按二进制进行编号，然后对各个编号按位倒置并按此重新排序。例如，对于一个8点变换，

001倒置以后变成 100

010 → 010

011 → 110

100 → 001

101 → 101

110 → 011

111 → 111

倒置后的编号为{0,4,2,6,1,5,3,7}。

然后将这n个点列作为输入传送到蝶形结网络中，注意将因子WNk{displaystyle W_{N}^{k}}逐层加入到蝶形网络中。

算法复杂度

由于按蝶形结网络计算n点变换要进行log n层计算，每层计算n个点的变换，故算法的时间复杂度为O(nlog⁡n){displaystyle mathrm {O} (nlog n)}。

其他算法

在Cooley-Tukey算法之外还有其他DFT的快速演算法。对于长度N=N1N2{displaystyle N=N_{1}N_{2}}且N1{displaystyle N_{1}}与N2{displaystyle 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看成是zN−1{displaystyle z^{N}-1}，并把它分解成zM−1{displaystyle z^{M-1}}与z2M+azM+1{displaystyle z^{2M}+az^{M}+1}的形式。

另一个从多项式观点的快速傅立叶变换法是Winograd算法。此演算法把zN−1{displaystyle z^{N}-1}分解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。这样只需要很少的乘法量（如果有需要的话），所以winograd是可以得到最少乘法量的快速傅立叶演算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd演算法让DFT可以用2k{displaystyle 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的结果会满足对称性：

XN−k{displaystyle mathbf {X} _{N-k}}=Xk∗{displaystyle mathbf {X} _{k}^{*}}

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

在Matlab上用一次N點FFT計算兩個N點實數信號x,y的FFT:

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

一度人们认为，用离散哈特利转换（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）转换


• 数论转换

複雜度以及運算量的極限

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

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

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

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

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

快速演算法查表:

當輸入信號長度為N時:

N = 1~60

N < 1000

N > 1000

参阅

• 离散傅里叶变换


• 并行快速傅里叶变换

参考资料

延伸阅读