快速傅立葉變換

  「FFT」重新導向至此。關於其他用法，請見「FFT (消歧義)」。

一個五項餘弦級數的時域訊號。頻率為 $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）或其反轉換的方法。傅立葉分析將訊號從原始域（通常是時間或空間）轉換到頻域的表示或者逆過來轉換。FFT會通過把DFT矩陣分解為稀疏（大多為零）因子之積來快速計算此類變換。 因此，它能夠將計算DFT的複雜度從只用DFT定義計算需要的 $O(n^{2})$，降低到 $O(n\log n)$，其中 $n$ 為數據大小。

快速傅立葉變換廣泛的應用於工程、科學和數學領域。這裏的基本思想在1965年才得到普及，但早在1805年就已推導出來。 1994年美國數學家吉爾伯特·斯特朗把FFT描述為「我們一生中最重要的數值演算法」，它還被IEEE科學與工程計算期刊列入20世紀十大演算法。

定義和速度

用FFT計算DFT會得到與直接用DFT定義計算相同的結果；最重要的區別是FFT更快。（由於捨入誤差的存在，許多FFT演算法還會比直接運用定義求值精確很多，後面會討論到這一點。）

令 x₀, ...., x_N-1 為複數。DFT由下式定義

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

直接按這個定義求值需要 O(N²) 次運算：X_k 共有 N 個輸出，每個輸出需要 N 項求和。直接使用DFT運算需使用N個複數乘法(4N 個實數乘法)與N-1個複數加法(2N-2個實數加法)，因此，計算使用DFT所有N點的值需要N²複數乘法與N²-N 個複數加法。FFT則是能夠在 O(N log N) 次操作計算出相同結果的任何方法。更準確的說，所有已知的FFT演算法都需要 O(N log N) 次運算（技術上O只標記上界），雖然還沒有已知的證據證明更低的複雜度是不可能的。

要說明FFT節省時間的方式，就得考慮複數相乘和相加的次數。直接計算DFT的值涉及到 N² 次複數相乘和 N(N−1) 次複數相加（可以通過削去瑣碎運算（如乘以1）來節省 O(N) 次運算）。眾所周知的基2庫利-圖基演算法，N 為2的冪，可以只用 (N/2)log₂(N) 次複數乘法（再次忽略乘以1的簡化）和 Nlog₂(N) 次加法就可以得到相同結果。在實際中，現代電腦通常的實際效能通常不受算術運算的速度和對複雜主體的分析主導，但是從 O(N²) 到 O(N log 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}}$

庫利-圖基演算法

主條目：庫利－圖基快速傅立葉變換演算法

庫利-圖基演算法是最常見的FFT演算法。這一方法以分治法為策略遞歸地將長度為$N=N_{1}N_{2}$的離散傅立葉變換分解為長度為$N_{1}$的$N_{2}$個較短序列的離散傅立葉變換，以及與$\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次單位根$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點變換。照這樣可繼續分解下去。這就是庫利-圖基快速傅立葉變換演算法的基本原理。根據主定理不難分析出此時演算法的時間複雜度為$\mathrm {O} (N\log N)$

演算法實現

• 蝶形結網絡和位反轉（Bit Reversal）：

首先將$n=2^{N}$個輸入點列按二進制進行編號，然後對各個編號按位元倒置並按此重新排序。例如，對於一個8點變換，

001倒置以後變成 100

000 → 000

001 → 100

010 → 010

011 → 110

100 → 001

101 → 101

110 → 011

111 → 111

倒置後的編號為{0,4,2,6,1,5,3,7}。

然後將這n個點列作為輸入傳送到蝶形結網絡中，注意將因子$W_{N}^{k}$逐層加入到蝶形網絡中。

演算法複雜度

由於按蝶形結網絡計算n點變換要進行log n層計算，每層計算n個點的變換，故演算法的時間複雜度為$\mathrm {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來表示（實數項為原本純實數資料的偶數項，虛數項則為奇數項）。

在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）轉換
• 數論轉換

複雜度以及運算量的極限

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

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

參閱

• 離散傅立葉變換
• 並列快速傅立葉變換
• 快速數論變換

延伸閱讀

• Brenner, N.; Rader, C. A New Principle for Fast Fourier Transformation. IEEE Acoustics, Speech & Signal Processing. 1976, 24 (3): 264–266. doi:10.1109/TASSP.1976.1162805. 
• Brigham, E. O. The Fast Fourier Transform. New York: Prentice-Hall. 2002. 
• Cormen, Thomas H.; Leiserson, Charles E.; Rivest, Ronald L.; Stein, Clifford. 30. (Polynomials and the FFT). Introduction to Algorithms 2. MIT Press and McGraw-Hill. 2001. ISBN 0-262-03293-7. 
• Duhamel, Pierre. Algorithms meeting the lower bounds on the multiplicative complexity of length-2ⁿ DFTs and their connection with practical algorithms. IEEE Trans. Acoust. Speech. Sig. Proc. 1990, 38 (9): 1504–151. doi:10.1109/29.60070. 
• Duhamel, P.; Vetterli, M. Fast Fourier transforms: a tutorial review and a state of the art. Signal Processing. 1990, 19: 259–299. doi:10.1016/0165-1684(90)90158-U. 
• Edelman, A.; McCorquodale, P.; Toledo, S. The Future Fast Fourier Transform?. SIAM J. Sci. Computing. 1999, 20: 1094–1114. doi:10.1137/S1064827597316266. 
• D. F. Elliott, & K. R. Rao, 1982, Fast transforms: Algorithms, analyses, applications. New York: Academic Press.
• Funda Ergün, 1995, , Proc. 27th ACM Symposium on the Theory of Computing: 407–416.
• Frigo, M.; Johnson, S. G. The Design and Implementation of FFTW3 (PDF). Proceedings of the IEEE. 2005, 93: 216–231 [2008-06-22]. doi:10.1109/jproc.2004.840301. （原始內容存檔 (PDF)於2019-07-16）. 
• H. Guo and C. S. Burrus, 1996, , Proc. SPIE Intl. Soc. Opt. Eng. 2825: 250–259.
• H. Guo, G. A. Sitton, C. S. Burrus, 1994, , Proc. IEEE Conf. Acoust. Speech and Sig. Processing (ICASSP) 3: 445–448.
• Steve Haynal and Heidi Haynal, "Generating and Searching Families of FFT Algorithms", Journal on Satisfiability, Boolean Modeling and Computation vol. 7, pp. 145–187 (2011).
• Heideman, Michael T.; Burrus, C. Sidney. On the number of multiplications necessary to compute a length-2ⁿ DFT. IEEE Trans. Acoust. Speech. Sig. Proc. 1986, 34 (1): 91–95. doi:10.1109/TASSP.1986.1164785. 
• Johnson, S. G.; Frigo, M. A modified split-radix FFT with fewer arithmetic operations (PDF). IEEE Trans. Signal Processing. 2007, 55 (1): 111–119 [2008-06-22]. doi:10.1109/tsp.2006.882087. （原始內容存檔 (PDF)於2021-02-25）. 
• 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.
• Kent, Ray D. and Read, Charles (2002). Acoustic Analysis of Speech. ISBN 978-0-7693-0112-9. Cites Strang, G. (1994)/May–June). Wavelets. American Scientist, 82, 250–255.
• Morgenstern, Jacques. Note on a lower bound of the linear complexity of the fast Fourier transform. J. ACM. 1973, 20 (2): 305–306. doi:10.1145/321752.321761. 
• Mohlenkamp, M. J. A fast transform for spherical harmonics (PDF). J. Fourier Anal. Appl. 1999, 5 (2–3): 159–184. doi:10.1007/BF01261607. （原始內容 (PDF)存檔於2007年2月6日）. 
• Nussbaumer, H. J. Digital filtering using polynomial transforms. Electronics Lett. 1977, 13 (13): 386–387. doi:10.1049/el:19770280. 
• 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).
• Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P., Chapter 12. Fast Fourier Transform, Numerical Recipes: The Art of Scientific Computing 3, New York: Cambridge University Press, 2007 [2015-12-11], ISBN 978-0-521-88068-8, （原始內容存檔於2011-08-11） 
• Rokhlin, Vladimir; Tygert, Mark. Fast algorithms for spherical harmonic expansions. SIAM J. Sci. Computing. 2006, 27 (6): 1903–1928. doi:10.1137/050623073. 
• Schatzman, James C. Accuracy of the discrete Fourier transform and the fast Fourier transform. SIAM J. Sci. Comput. 1996, 17: 1150–1166. doi:10.1137/s1064827593247023. 
• Shentov, O. V.; Mitra, S. K.; Heute, U.; Hossen, A. N. Subband DFT. I. Definition, interpretations and extensions. Signal Processing. 1995, 41 (3): 261–277. doi:10.1016/0165-1684(94)00103-7. 
• Sorensen, H. V.; Jones, D. L.; Heideman, M. T.; Burrus, C. S. Real-valued fast Fourier transform algorithms. IEEE Trans. Acoust. Speech Sig. Processing. 1987, 35 (35): 849–863. doi:10.1109/TASSP.1987.1165220.  See also Sorensen, H.; Jones, D.; Heideman, M.; Burrus, C. Corrections to "Real-valued fast Fourier transform algorithms". IEEE Transactions on Acoustics, Speech, and Signal Processing. 1987, 35 (9): 1353–1353. doi:10.1109/TASSP.1987.1165284. 
• Welch, Peter D. A fixed-point fast Fourier transform error analysis. IEEE Trans. Audio Electroacoustics. 1969, 17 (2): 151–157. doi:10.1109/TAU.1969.1162035. 
• Winograd, S. On computing the discrete Fourier transform. Math. Computation. 1978, 32 (141): 175–199. JSTOR 2006266. doi:10.1090/S0025-5718-1978-0468306-4.