本頁使用了標題或全文手工轉換

快速傅立葉變換

維基百科,自由的百科全書
前往: 導覽搜尋

快速傅立葉變換英語Fast Fourier Transform, FFT),是離散傅立葉變換的快速演算法,也可用於計算離散傅立葉變換的逆變換。快速傅立葉變換有廣泛的應用,如數位訊號處理、計算大整數乘法、求解偏微分方程等等。本條目只描述各種快速演算法。

對於複數序列x_{0},\ x_{1},\ ...,\ x_{n-1},離散傅立葉變換公式為:


 y_j = \sum_{k=0}^{n-1} e^{-{2\pi\imath \over n} j k}x_k \qquad j = 0,1,\dots,n-1.

直接變換的計算複雜度是\Omicron(n^2)(參見大O符號)。快速傅立葉變換可以計算出與直接計算相同的結果,但只需要\Omicron (n \log n)的計算複雜度。通常,快速演算法要求n能被因數分解,但不是所有的快速傅立葉變換都要求n合數,對於所有的整數n,都存在複雜度為\Omicron (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}^TM_0個相異的非平凡值(a_m\ne\pm2^k,a_m\ne\pm2^ka_n where  m\ne 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_1X[1]+b_2X[2]+\cdots+b_nX[N]

Step 2:Z[1]=a_1Z_a,Z[2]=a_2Z_a,\cdots,Z[N]=a_mZ_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_1P_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_2B_1+N_1B_2+3D_1+2D_2}

庫利-圖基演算法[編輯]

庫利-圖基演算法是最常見的FFT演算法。這一方法以分治法為策略遞歸地將長度為N=N_1 N_2離散傅立葉變換分解為長度為N_1N_2個較短序列的離散傅立葉變換,以及與\Omicron (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演算法聯合使用。

目前最常用的快速傅立葉變換是庫利-圖基(Cooley-Tukey)演算法。該演算法利用分治法,將一個常數n的離散傅立葉變換遞歸地分解為兩個常數分別為n_1n_2的變換,保證n=n_1 n_2,從而簡化了原來的離散傅立葉變換。

設計思想[編輯]

下面,我們用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. mN的約數,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點變換。照這樣可繼續分解下去。這就是庫利-圖基快速傅立葉變換演算法的基本原理。根據主定理不難分析出此時演算法的時間複雜度為\Omicron (N\log N)


演算法實現[編輯]

  • 蝶形結網路和位反轉(Bit Reversal):
    • 首先將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個點列作為輸入傳送到蝶形結網路中,注意將因子W_{N}^{k}逐層加入到蝶形網路中。

演算法複雜度[編輯]

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

其他演算法[編輯]

Cooley-Tukey演算法之外還有其他DFT的快速演算法。對於長度N = N_1N_2N_1N_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可用下列的演算法代替:

複雜度以及運算量的極限[編輯]

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

在1978年,Winograd率先導出一個較嚴謹的FFT所需乘法量的下限:\Theta (N)。當N=2^k時,DFT只需要4N-2\log_{2}^2N-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)。

參考資料[編輯]

參閱[編輯]