檢視快速傅立葉變換的原始碼

'''快速傅立葉變換'''（英語：'''Fast Fourier Transform , FFT'''）， 是快速計算序列的離散傅立葉變換（DFT）抑是講反轉換的方法。傅立葉分析共訊號對原始域（通常是時間抑是空間）轉換到頻域的表示抑是倒反的轉換。FFT 會通過共 DFT 矩陣分解為疏櫳（大多為零）因為積來快速計算這類變換。所以，伊會當將計算 DFT 的複雜度對干焦用 DFT 定義計算需要的 $ O ( n ^ { 二 } ) $，降低到 $ O ( n \ log n ) $，其中 $ n $ 為資料大細。

快速傅立葉變換廣泛的應用佇工程、科學佮數學領域。遮的基本思想佇一九六五年才得著普遍，毋過早佇一八零五年就已經推出來。一九九四年美國數學家吉爾伯特 ・ 斯特朗共 FFT 欲描述講「咱一生中上重要的數值演算法」，伊閣予人 IEEE 科學和工程計算期刊列入二十世紀十大演算法。

==定義佮速度==

用 FFT 計算 DFT 得著佮直接用 DFT 定義計算相仝的結果；上重要的區別是 FFT 較緊咧。（因為重耽入精差的存在，真濟 FFT 演算法猶閣會比直接運用定義求值精確真濟，後壁會討論著這點仔。）

令 _ x _ 零 , . . . . , _ x _ N 交一為複數。DFT 由下式定義


: $ X _ { k }=\ sum _ { n=零 } ^ { N 影一 } x _ { n } e ^ {-{ i 二 \ pi k { \ frac { n } { N } } } } \ qquad k=零 , \ dots , N 影一 . $

直接揤這个定義求值需要 O ( _ N _ 二 ) 次運算：_ X _ k 共有 _ N _ 喝一下，每一个輸出需要 _ N _ 項求和。直接使用 DFT 運算需要使用 N 個複數乘法 ( 四 N 個實數乘法 ) 佮 N 學一个複數加法 ( 二 N 學二个實數加法 )，所以，算使用 DFT 所有 N 點的值需要 _ N _ 二複數乘法佮 _ N _ 二-N 個複數加法。FFT 著會當佇 O ( _ N _ log _ N _ ) 操作計算出仝款結果的任何方法。閣較準確實講，所有已經知的 FFT 演算法攏需要 O ( _ N _ log _ N _ ) 次運算（技術上 O 干焦標記牢界）， 雖然猶無已知的證據證明閣較低的複雜度是無可能的。

愛說明 FFT 愛節省時間的方式，就愛考慮複數相乘佮相加的次數。直接計算 DFT 的值牽連著 _ N _ 二次複數相乘和 _ N _ ( _ N _ − 一 ) 次複數相加（會使通過削去工夫碎運算（若乘以一）來節省 O ( _ N _ ) 次運算）。 眾所周知的基二庫利-圖基演算法，_ N _ 為二的冪，會使干焦用 ( _ N _ / 二 ) log 二 ( _ N _ ) 次複數乘法（閣再忽略乘以一的簡化）和 _ N _ log 二 ( _ N _ ) 次加法就會用得著仝款結果。佇實際上內底，現代電腦通常的實際效能通常無受著算術運算的速度佮對複雜主體的分析主導，猶毋過對 O ( _ N _ 二 ) 到 O ( _ N _ log _ N _ ) 的母體改進猶原會當體出來。

==一般的簡單化理論==

準講一个 M \ * N 型矩陣'''S'''會當分解做列向量佮行向量相乘：

$ $
\ mathbf { S }={ \ begin { bmatrix } a _ { 一 } \ \ a _ { 二 } \ \ \ vdots \ \ a _ { m } \ end { bmatrix } } { \ begin { bmatrix } b _ { 一 } & b _ { 二 } & \ cdots & b _ { n } \ end { bmatrix } }
$ $

若是 $ { \ begin { bmatrix } a _ { 一 } & a _ { 二 } & \ cdots & a _ { m } \ end { bmatrix } } ^ { T } $ 有 $ M _ { 零 } $ 一个欲相異的非平凡值（$ a _ { m } \ neq \ pm 二 ^ { k } , a _ { m } \ neq \ pm 二 ^ { k } a _ { n } $ where $ m \ neq n $）

$ { \ begin { bmatrix } b _ { 一 } & b _ { 二 } & \ cdots & b _ { n } \ end { bmatrix } } $ 有 $ N _ { 零 } $ 個欲相異的非平凡值則'''S'''共需要 $ M _ { 零 } * N _ { 零 } $ 個乘法。

$ $
{ \ begin { bmatrix } Z [一] \ \ Z [二] \ \ \ vdots \ \ Z [N] \ end { bmatrix } }=\ mathbf { S } { \ begin { bmatrix } X [一] \ \ X [二] \ \ \ vdots \ \ X [N] \ end { bmatrix } }={ \ begin { bmatrix } a _ { 一 } \ \ a _ { 二 } \ \ \ vdots \ \ a _ { m } \ end { bmatrix } } { \ begin { bmatrix } b _ { 一 } & b _ { 二 } & \ cdots & b _ { n } \ end { bmatrix } } { \ begin { bmatrix } X [一] \ \ X [二] \ \ \ vdots \ \ X [N] \ end { bmatrix } }
$ $

Step 一：$ Z _ { a }=b _ { 一 } X [一] + b _ { 二 } X [二] + \ cdots + b _ { n } X [N] $

Step 二：$ Z [一]=a _ { 一 } Z _ { a } , Z [二]=a _ { 二 } Z _ { a } , \ cdots , Z [N]=a _ { m } Z _ { a } $

簡化理論的變型：

$ $
\ mathbf { S }={ \ begin { bmatrix } a _ { 一 } \ \ a _ { 二 } \ \ \ vdots \ \ a _ { m } \ end { bmatrix } } { \ begin { bmatrix } b _ { 一 } & b _ { 二 } & \ cdots & b _ { n } \ end { bmatrix } } + \ mathbf { S } _ { 一 }
$ $

$ S _ { 一 } $ 嘛是一个 M \ * N 矩陣。

若是 $ S _ { 一 } $ 有 $ P _ { 一 } $ 個值無等於零，著 $ \ mathbf { S } $ 的乘法量上限為 $ M _ { 零 } + N _ { 零 } + P _ { 一 } $。

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

準講 $ N=P _ { 一 } \ times P _ { 二 } \ times \ cdots \ times P _ { k } $，其中 $ P _ { 一 } , P _ { 二 } , \ cdots , P _ { k } $ 互相落氣

$ \ mathbf { P _ { k } } $ 點 DFT 的乘法量為 $ \ mathbf { B _ { k } } $，著 $ \ mathbf { N } $ 點 DFT 的乘法量為：


: $ { \ frac { N } { P _ { 一 } } } B _ { 一 } + { \ frac { N } { P _ { 二 } } } B _ { 二 } + \ cdots \ cdots + { \ frac { N } { P _ { k } } } B _ { k } $

準講 $ \ mathbf { N=P ^ { c } } $，P 是一个質數。

若是 $ \ mathbf { N _ { 一 }=P ^ { a } } $ 點的 DFT 需要的乘法量為著 $ \ mathbf { B _ { 一 } } $

而且 $ n _ { 一 } \ times n _ { 二 } $ 當中（$ n _ { 一 }=零 , 一 , \ cdots , N _ { 一 } 影一 , \ quad n _ { 二 }=零 , 一 , \ cdots , N _ { 二 } 影一 $）

有 $ D _ { 一 } $ 個價值不為 $ { \ frac { N } { 十二 } } $ 佮 $ { \ frac { N } { 八 } } $ 的倍數，

有 $ D _ { 二 } $ 個值為 $ { \ frac { N } { 十二 } } $ 佮 $ { \ frac { N } { 八 } } $ 的倍數，但是不為講 $ { \ frac { N } { 四 } } $ 的倍數，

著 N 點 DFT 的乘法量為：


: $ \ mathbf { N _ { 二 } B _ { 一 } + N _ { 一 } B _ { 二 } + 三 D _ { 一 } + 二 D _ { 二 } } $

==庫利-圖基演算法==

庫利-圖基演算法是上捷看著的 FFT 演算法。這一方法以分治法為策略遞迴地會長度為 $ N=N _ { 一 } N _ { 二 } $ 的離散傅立葉變換分解做長度為 $ N _ { 一 } $ 的 $ N _ { 二 } $ 較短序列的離散傅立葉變換，以及佮 $ \ mathrm { O } ( N ) $ 個旋轉因為的複數乘法。

這種方法佮 FFT 基本想講佇一九六五年 J . W . Cooley 和 J . W . Tukey 合作發表 _ An algorithm for the machine calculation of complex Fourier series _ 了開始共人知影講。但是後來發現講，實際上這兩位作者只是重新發明矣高斯一八空五年就已經提出的演算法（此演算法佇歷史上幾若改以各種的形式被再一改提出）。

庫利-圖基演算法上出名的應用，是將序列長做 _ N _ 的 DFT 分割為兩个長為 _ N / 二 _ 的子序列的 DFT，所以這應用干焦適用佇序列長度做二的冪的 DFT 計算，即基二-FFT。實際上，親像高斯和 Cooley 佮 Tukey 攏指出的彼款的，Cooley-Tukey 演算法嘛會當產列長度 _ N _ 為任意因數分解形式的 DFT，即混合基 FFT，而且閣會使應用佇其他的諸如分裂基 FFT 等變種。就算講 Cooley-Tukey 演算法的基本思路是用遞迴的方法進行計算，大多數傳統的演算法實現攏將顯式的遞迴演算法治解寫做非遞迴的形式。另外咧，因為乎 Cooley-Tukey 演算法是將 DFT 分解為較細長度的濟个 DFT，因此伊會當仝任一種其他的 DFT 演算法聯合使用。

===設計思想===

下跤，阮用'''N 次單位根'''$ W _ { N } $ 來表示 $ e ^ {-j { \ frac { 二 \ pi } { N } } } $。

$ W _ { N } $ 的性質：

一 .'''週期性'''，$ W _ { N } $ 有禮拜 $ N $，即 $ W _ { N } ^ { k + N }=W _ { N } ^ { k } $
二 .'''對稱性'''：$ W _ { N } ^ { k + { \ frac { N } { 二 } } }=-W _ { N } ^ { k } $。
三 . 若是 $ m $ 是 $ N $ 的因數，$ W _ { N } ^ { mkn }=W _ { \ frac { N } { m } } ^ { kn } $

為著簡單起見，咱下跤設待變換序列長度 $ n=二 ^ { r } $。根據頂頭單位根的對稱性，求級數 $ y _ { k }=\ sum _ { n=零 } ^ { N 影一 } W _ { N } ^ { kn } x _ { n } $ 時，會當將求佮區間來分做兩部份：

$ $
{ \ begin { matrix } y _ { k }=\ sum _ { n=二 t } W _ { N } ^ { kn } x _ { n } + \ sum _ { n=二 t + 一 } W _ { N } ^ { kn } x _ { n } \ \=\ sum _ { t } W _ { \ frac { N } { 二 } } ^ { kt } x _ { 二 t } + W _ { N } ^ { k } \ sum _ { t } W _ { \ frac { N } { 二 } } ^ { kt } x _ { 二 t + 一 } \ \=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 \ } _ { 零 } ^ { N 影一 } $ 奇數號佮偶數號序列 N / 二點變換。由此式只會當計算出 $ y _ { k } $ 的前 N / 兩个點，對於後 N / 兩个點，注意 $ F _ { odd } ( k ) $ 和 $ F _ { even } ( k ) $ 攏是禮拜為 N / 二的函式，由單位根的對稱性，所以變換公式的啦：

* $ y _ { k + { \ frac { N } { 二 } } }=F _ { even } ( k )-W _ { N } ^ { k } F _ { odd } ( k ) $
* $ y _ { k }=F _ { even } ( k ) + W _ { N } ^ { k } F _ { odd } ( k ) $。

按呢乎，一个 N 點變換就分解講兩个 N / 二點變換。照按呢會當繼續分解落去。這就是'''庫利-圖基快速傅立葉變換'''演算法的基本原理。根據主定理袂歹分析出現佇這个時陣佇算法的時間的複雜度做 $ \ mathrm { O } ( N \ log N ) $

===演算法實現===

* 蝶仔形結網路佮位反轉（Bit Reversal）：


: 首先將 $ n=二 ^ { N } $ 輸入點列揤二進位進行編號，然後對逐个編號照位元倒置閣按呢重新鋪排。比如講，對一个八點變換，


: 一倒置了後變成百外


: 零 → 零


: 一 → 一百


: 十 → 十


: 十一 → 一百十一


: 一百 → 一


: 一百空一 → 一百空一


: 一百十一 → 十一


: 一百十一 → 一百十一


: 倒反的時陣的編號做 {零 , 四 , 二 , 六 , 一 , 五 , 三 , 七 }。


: 才閣將這 n 一个點列做輸入傳送到蝶形結網路內底，注意共因為 $ W _ { N } ^ { k } $ 每一層加入到蝶形網路中。

===演算法複雜度===

因為按蝶仔形結網路算 n 點變換欲做 log _ n _ 層計算，逐層計算 n 個點的變換，故演算法的時間複雜度做 $ \ mathrm { O } ( n \ log n ) $。

==其他演算法==

佇咧 Cooley-Tukey 算法之外閣有其他 DFT 的快速演算法。對著長度 $ N=N _ { 一 } N _ { 二 } $ 而且 $ N _ { 一 } $ 佮 $ N _ { 二 } $ 互質的序列，會當採用是中國賰定理的互質因為子的演算法共 N 長序列的 DFT 分解為兩个囝序列的 DFT。佮 Cooley-Tukey 演算法無仝的是，互質因為演算法無需要轉踅因為。

Rader-Brenner 演算法是類似 Cooley-Tukey 演算法，但是用的旋轉因為攏是純虛數，以增加法運算和降低了數值穩定性做代價減少了乘法運算。這方法了後去予 split-radix variant of Cooley-Tukey 所取代，佮 Rader-Brenner 演算法相比，有平濟的乘法量，煞有較少的加法量，而且無犧牲數值的準確性。

Bruun 以及 QFT 演算法是不斷的共 DFT 分做真濟較細的 DFT 運算。（Rader-Brenner 以及 QFT 演算法是為著二的指數所設計的演算法，但是猶原會當適用佇咧會當分解的整數頂頭。Bruun 演算法則會使運用佇會使被分做偶數個運算的數位）。 尤其是 Bruun 演算法，共 FFT 看做是 $ z ^ { N } 影一 $，閣共分解做 $ z ^ { M 影一 } $ 佮 $ z ^ { 二 M } + az ^ { M } + 一 $ 彼个形體。

另外一个對濟項式觀點的快速傅立葉變換法是 Winograd 演算法。此演算法共 $ z ^ { N } 影一 $ 分解做 cyclotomic 多項式，這多項式的係數通常為一，零，影一。按呢只要有足少的乘法量（若是有需要）， 所以乎 winograd 是會當得著上少乘法量的快速傅立葉演算法，對著較細項的數位，會當揣出有效率的算方式。閣較精確咧講，winograd 演算法予 DFT 會用得 $ 二 ^ { k } $ 點的 DFT 來簡化，但是減少乘法量的同時，嘛增加足濟的加法量。Winograd 嘛會當利用賰的值定理來簡化 DFT。

Rader 演算法提出了利用點數為 N（N 為質數）的 DFT 進行長度為 N 鋪一的迴旋拗積來表示阮原本的 DFT，按呢就會當利用拗積用一對基本的 FFT 來計算 DFT。另外一个 prime-size 的 FFT 演算法為 chirp-Z 演算法。此法嘛是將 DFT 用拗積來表示，此法佮 Rader 演算法相比，會當運用轉踅，其轉換的基礎為 Z 轉換（Rabiner et al . , 一千九百六十九）。

==實數抑是對稱資料專用的演算法==

佇真濟的運用當中，欲來進行 DFT 的資料是純實數，按呢經過 DFT 結果會滿足對稱性：


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

佇咧有一寡演算法是專門為這種情形設計的（e . g . Sorensen , 一千九百八十七）。 另外一寡利用舊的演算法（e . g . Cooley-Tukey）， 刪去一寡無必要的去做算步，按呢省有記持的使用，原仔省時間。另外一方面，嘛會當共一个偶數長度而且純實數的 DFT，用長度為原本一半的複數型態 DFT 來表示（實數項為原本純實數資料的偶數項，虛數項為奇數項）。

'''佇咧 Matlab 上用一改 N 點 FFT 計算兩个 N 點實數訊號 x , y 的 FFT :'''

一度人認為講，用離散哈特利換（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）轉換
* 數論轉換

==複雜度佮運算量的極限==

長久以來，人對求出快速傅立葉變換的複雜度下限以及需要偌濟的運算量感覺誠有興趣，實際上嘛閣有真濟問題需要解決。就算是用較簡單的情形，即 $ 二 ^ { k } $ 點的 DFT，嘛猶無法度嚴謹的證明出 FFT 上無需要 $ \ Omega ( N \ log N ) $（不比 $ N \ log N $ 細）的運算量，目前嘛無發現複雜度閣較低的演算法。通常數學運算量的濟少就是運算效率好䆀上主要的因素，但是佇現實中，有真濟因素嘛會有真大的影響，若緊取記持體以及 CPU 攏有足大的影響。

佇一九七八年，Winograd 率先導出一个較嚴謹的 FFT 所需乘法量的下限：$ \ Theta ( N ) $。當 $ N=二 ^ { k } $ 時，DFT 只需要 $ 四 N 鋪二 \ log _ { 二 } ^ { 二 } N 鋪二 \ log _ { 二 } N 扳四 $ 次無理實數的乘法即會當計算出來。閣較詳盡，而且嘛會當近這下限的演算法嘛一一予人提出（Heideman & Burrus , 一千九百八十六 ; Duhamel , 一千九百九十）。 足可惜的是，遮的演算法，攏需要足大量的加法計算，目前的硬體無法度克服這个問題。

對這个所需要加法量的數目，雖然咱會當佇某一寡受限制的假影，捒甲其下限，但目前並無一个精確的下限去予人推導出來。一九七三年，Morgenstern 佇乘法常數趨近誠大的情形之下（著大部份的 FFT 演算法為真，但是毋是全部）推導出加法量的下限：$ \ Omega \ left ( N \ log N \ right ) $。Pan（一千九百八十六）佇咧假使 FFT 演算法的無仝步的情形有其極限下證明出加法量的下限 $ \ Omega ( NlogN ) $，但是一般來講，此假使相當的不明確。長度為 $ N=二 ^ { k } $ 的情形下，佇某一寡假影，Papadimitriou（一千九百七十九）提出使用 Cooley-Tukey 演算法所需要的複數加法量 $ N \ log _ { 二 } N $ 上少的是。到這當時為止，佇咧長度為 $ N=二 ^ { k } $ 狀況，猶無任何 FFT 的演算法會當予複數的加法比 $ N \ log _ { 二 } N $ 猶閣少。

猶閣有一个問題是按怎共乘法量佮加法量的總和最小化，有當時仔號做 " 演算複雜度 "（佇遮考慮的是實際的運算量，毋是漸漸複雜度）。 按呢仝款，無一个頂真下限去予證明出來。對一九六八年開始，$ N=二 ^ { k } $ 點 DFT 來講，split-radix FFT 演算法需要上少的運算量，佇咧 $ N > 一 $ 的情形下，其需要 $ 四 N \ log _ { 二 } N ma六 N + 八 $ 個乘法運算以及加法運算。最近有人導出閣較低的運算量：$ { \ frac { 三十四 } { 九 } } N \ log _ { 二 } N $。（Johnson and Frigo , 兩千空七 ; Lundy and Van Buskirk , 兩千空七）

大多數試欲降低或者是證明 FFT 複雜度下限的人攏共焦點囥佇咧複數資料輸入的情況，因為上簡單的情形。猶毋過，複數資料輸入的 FFT 演算法，佮實數資料輸入的 FFT 演算法，離散餘弦轉換（DCT）， 離散哈特列轉換（DHT）， 猶閣有其他的演算法，攏有足大的關連性。故任何一个演算法，佇複雜度上有任何的改善，其他的演算法複雜度嘛會隨得著改善（Duhamel & Vetterli , 一千九百九十）。

==快速演算法查表==

做彼个輸入訊號長度做 N 時 :

N=一 ~ 六十

N < 一千

N > 一千

==參閱==

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

==參考資料==

==延伸閱讀==

[[分類: 待校正]]