快速傅立葉變換

快速計算序列的離散傅立葉變換或其逆變換的方法

快速傅立葉變換(英語:Fast Fourier Transform, FFT),是快速計算序列的離散傅立葉變換(DFT)或其反轉換的方法[1]傅立葉分析將訊號從原始域(通常是時間或空間)轉換到頻域的表示或者逆過來轉換。FFT會通過把DFT矩陣分解稀疏(大多為零)因子之積來快速計算此類變換。[2] 因此,它能夠將計算DFT的複雜度從只用DFT定義計算需要的 ,降低到 ,其中 為資料大小。

一個五項餘弦級數的時域訊號。頻率為 的倍數。振幅為 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),代價是速度慢了上千倍。

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

定義和速度

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

x0, ...., xN-1複數。DFT由下式定義

 

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

要說明FFT節省時間的方式,就得考慮複數相乘和相加的次數。直接計算DFT的值涉及到 N2 次複數相乘和 N(N−1) 次複數相加(可以通過削去瑣碎運算(如乘以1)來節省 O(N) 次運算)。眾所周知的基2庫利-圖基演算法N 為2的冪,可以只用 (N/2)log2(N) 次複數乘法(再次忽略乘以1的簡化)和 Nlog2(N) 次加法就可以得到相同結果。在實際中,現代電腦通常的實際效能通常不受算術運算的速度和對複雜主體的分析主導[7],但是從 O(N2) 到 O(N log N) 的母體改進仍然能夠體現出來。

一般的簡化理論

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

 

  個相異的非平凡值(  where  ) 

  個相異的非平凡值

S共需要 個乘法。

 

Step 1: 

Step 2: 

簡化理論的變型:

 

 也是一個M*N的矩陣。

  個值不等於0,則 的乘法量上限為 

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

假設 ,其中 彼此互質

 點DFT的乘法量為 ,則 點DFT的乘法量為:

 

假設 ,P是一個質數。

 點的DFT需要的乘法量為 

 當中( 

 個值不為  的倍數,

 個值為  的倍數,但不為 的倍數,

則N點DFT的乘法量為:

 

庫利-圖基演算法

庫利-圖基演算法是最常見的FFT演算法。這一方法以分治法為策略遞迴地將長度為 離散傅立葉變換分解為長度為  個較短序列的離散傅立葉變換,以及與 旋轉因子的複數乘法。

這種方法以及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次單位根 來表示 

 的性質:

  1. 週期性 具有週期 ,即 
  2. 對稱性 
  3.   的因數, 

為了簡單起見,我們下面設待變換序列長度 。根據上面單位根的對稱性,求級數 時,可以將求和區間分為兩部分:

 

  是兩個分別關於序列 奇數號和偶數號序列N/2點變換。由此式只需計算出 的前N/2個點,對於後N/2個點可以通過複指數函式的對稱性快速求解。即,注意  都是週期為N/2的函式,由單位根的對稱性,有以下變換公式:

  •  
  •  

這樣,一個N點變換就分解成了兩個N/2點變換。照這樣可繼續分解下去。這就是庫利-圖基快速傅立葉變換演算法的基本原理。根據主定理不難分析出此時演算法的時間複雜度為 

演算法實現

  • 蝶形結網路和位反轉(Bit Reversal):
首先將 個輸入點列按二進制進行編號,然後對各個編號按位元倒置並按此重新排序。例如,對於一個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個點列作為輸入傳送到蝶形結網路中,注意將因子 逐層加入到蝶形網路中。

演算法複雜度

由於按蝶形結網路計算n點變換要進行log n層計算,每層計算n個點的變換,故演算法的時間複雜度為 

其他演算法

Cooley-Tukey演算法之外還有其他DFT的快速演算法。對於長度   互質的序列,可以採用基於中國剩餘定理互質因子演算法將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看成是 ,並把它分解成  的形式。

另一個從多項式觀點的快速傅立葉變換法是Winograd演算法。此演算法把 分解成cyclotomic多項式,而這些多項式的係數通常為1,0,-1。這樣只需要很少的乘法量(如果有需要的話),所以winograd是可以得到最少乘法量的快速傅立葉演算法,對於較小的數位,可以找出有效率的算方式。更精確地說,winograd演算法讓DFT可以用 點的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的結果會滿足對稱性:

 = 

而有一些演算法是專門為這種情形設計的(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可用下列的演算法代替:

複雜度以及運算量的極限

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

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

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

還有一個問題是如何把乘法量與加法量的總和最小化,有時候稱作"演算複雜度"(在這裡考慮的是實際的運算量,而不是漸近複雜度)。同樣的,沒有一個嚴謹下限被證明出來。從1968年開始, 點DFT而言,split-radix FFT演算法需要最少的運算量,在 的情形下,其需要 個乘法運算以及加法運算。最近有人導出更低的運算量: 。(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

參閱

參考資料

  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

延伸閱讀