Difference between revisions of "快速傅立叶变换"

From SoundDB
Jump to: navigation, search
(Created page with "{{noteTA |T=zh-hans:快速傅里叶变换; zh-hant:快速傅立葉轉換; |G1=Communication|1=zh:傅里叶; zh-hans:傅里叶; zh-hant:傅立葉; }} {{傅里叶变换}} '''快速...")
 
(No difference)

Latest revision as of 01:24, 6 February 2011

Template:NoteTA Template:傅里叶变换 快速傅里叶变换（Template:Lang，Template:Lang），是离散傅里叶变换的快速算法，也可用于计算离散傅里叶变换的逆变换。快速傅里叶变换有广泛的应用，如数字信号处理、计算大整数乘法、求解偏微分方程等等。本条目只描述各种快速算法，对于离散傅里叶变换的性质和应用，请参见离散傅里叶变换。

对于复数序列<math>x_{0},\ x_{1},\ ...,\ x_{n-1}</math>，离散傅里叶变换公式为：

<math> 

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

</math>

直接变换的计算复杂度是<math>\mathcal{O}(n^2)</math>（参见大O符号）。快速傅里叶变换可以计算出与直接计算相同的结果，但只需要<math>\mathcal{O}(n \log n)</math>的计算复杂度。通常，快速算法要求n能被因数分解，但不是所有的快速傅里叶变换都要求n是合数，对于所有的整数n，都存在复杂度为<math>\mathcal{O}(n \log n)</math>的快速算法。

除了指数的符号相反、并多了一个1/n的因子，离散傅里叶变换的正变换与逆变换具有相同的形式。因此所有的离散傅里叶变换的快速算法同时适用于正逆变换。

一般的簡化理論

假設一個M*N Sub-rectangular matrix S可分解成列向量以及行向量相乘：

<math>\mathbf{S}=\begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_n\end{bmatrix}\begin{bmatrix} b_1 & b_2 & \cdots & b_n\end{bmatrix}</math>

若<math>\begin{bmatrix} a_1 & a_2 & \cdots & a_n\end{bmatrix}^T</math>有<math>M_0</math>個相異的non-trivial values(<math>a_m\ne\pm2^k,a_m\ne\pm2^ka_n </math> where <math> m\ne n</math>)

<math>\begin{bmatrix} b_1 & b_2 & \cdots & b_n\end{bmatrix}</math>有<math>N_0</math>個相異的non-trivial values

則S共需要<math>M_0+N_0</math>個乘法。

<math>\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_n\end{bmatrix}\begin{bmatrix} b_1 & b_2 & \cdots & b_n\end{bmatrix}\begin{bmatrix} X[1] \\ X[2] \\ \vdots \\ X[N]\end{bmatrix}</math>

Step 1：<math>Z_a=b_1X[1]+b_2X[2]+\cdots+b_nX[N]</math>

Step 2：<math>Z[1]=a_1Z_a,Z[2]=a_2Z_a,\cdots,Z[N]=a_nZ_a</math>

簡化理論的變型：

<math>\mathbf{S}=\begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_n\end{bmatrix}\begin{bmatrix} b_1 & b_2 & \cdots & b_n\end{bmatrix}+\mathbf{S}_1</math>

<math>S_1</math>也是一個M*N的矩陣。

若<math>S_1</math>有<math>P_1</math>個值不等於0，則<math>\mathbf{S}</math>的乘法量上限為<math>M_0+N_0+P_1</math>。

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

假設<math>N= P_1 \times P_2 \times \cdots \times P_k \qquad P_1,P_2, \cdots , P_k </math>彼此互質

<math>\mathbf{P_k}</math> 點DFT的乘法量為<math> \mathbf{B_k}</math>，則<math>\mathbf{N}</math>點DFT的乘法量為：

<math>\frac{N}{P_1}B_1+\frac{N}{P_2}B_2+\cdots\cdots+\frac{N}{P_k}B_k</math>

假設<math>\mathbf{N=P^c}</math>，P是一個質數。

若<math>\mathbf{N_1=P^a}</math>點的DFT需要的乘法量為<math>\mathbf{B_1}</math>

且<math>n_1\times n_2</math> 當中 (<math> n_1=0,1,\cdots,N_1-1, \quad n_2=0,1, \cdots , N_2-1</math>)

有<math>D_1</math>個值不為<math>\frac{N}{12}</math>及<math>\frac{N}{8}</math>的倍數，

有<math>D_2</math>個值為<math>\frac{N}{12}</math>及<math>\frac{N}{8}</math>的倍數，但不為<math>\frac{N}{4}</math>的倍數，

則N點DFT的乘法量為：

<math>\mathbf{N_2B_1+N_1B_2+3D_1+2D_2}</math>

Cooley-Tukey算法

Template:Main

Cooley-Tukey算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为<math>N=N_1 N_2</math>的DFT分解为长度分别为<math>N_1</math>和<math>N_2</math>的两个较短序列的DFT，以及与<math>\mathcal{O}(N)</math>个旋转因子的复数乘法。

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

Cooley-Tukey算法最有名的应用，是将序列长为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次单位根<math>W_{N}</math>来表示<math>e^{-j\frac{2\pi}{N}}</math>。

<math>W_{N}</math>的性质：

  1. 周期性，<math>W_{N}</math>具有周期N。
  2. 对称性：<math>W_{N}^{k+\frac{N}{2}}=-W_{N}^{k}</math>。
  3. <math>W_{N}^{ikn}=W_{\frac{N}{i}}^{kn}</math>

为了简单起见，我们下面设待变换序列长度<math>n=2^r</math>。 根据上面单位根的对称性，求级数<math>y_k=\sum_{n=0}^{N-1} W_{N}^{kn}x_n</math>时，可以将求和区间分为两部分：

<math>\begin{matrix}y_k=\sum_{n=2i} W_{N}^{kn} x_n + \sum_{n=2i+1} W_{N}^{kn}x_n\\= \sum_{i} W_{\frac{N}{2}}^{ki}x_{2i} + W_{N}^{k}\sum_{i} W_{\frac{N}{2}}^{ki}x_{2i+1}\\= F_{even}(k) + W_{N}^{k}F_{odd}(k)&&&&&&(i\in\mathbb{Z})\end{matrix}</math>

<math>F_{odd}(k)</math> 和 <math>F_{even}(k)</math>是两个分别关于序列<math>\left\{x_n\right\}_0^{N-1}</math>奇数号和偶数号序列N/2点变换。由此式只能计算出<math>y_k</math>的前N/2个点，对于后N/2个点，注意 <math>F_{odd}(k)</math> 和 <math>F_{even}(k)</math> 都是周期为N/2的函数，由单位根的对称性，于是有以下变换公式：

  • <math>y_{k+\frac{N}{2}} = F_{even}(k) - W_{N}^{k}F_{odd}(k)</math>
  • <math>y_k = F_{even}(k) + W_{N}^{k}F_{odd}(k)</math>。

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


其他算法

在Cooley-Tukey算法之外還有其他DFT的快速演算法。對於長度N = N1N2且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看成是<math>z^N-1</math>，並把它分解成<math>z^{M-1}</math> 與<math>z^{2M}+az^M+1</math> 的形式。

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

<math>\mathbf{X}_{N-k}</math>=<math>\mathbf{X}_k^*</math>

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

  • DCT.
  • DST.
  • DHT.
  • 直接做摺積(direct computing)
  • 分段式DFT摺積(sectioned DFT convolution)
  • Winograd 演算法
  • Walsh(Hadamard)轉換
  • 数论轉換

複雜度以及運算量的極限

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

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

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

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

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

参考资料

参阅

ca:Transformada Ràpida de Fourier cs:Rychlá Fourierova transformace da:Fast Fourier Transform de:Schnelle Fourier-Transformation en:Fast Fourier transform es:Transformada rápida de Fourier fa:تبدیل سریع فوریه fr:Transformée de Fourier rapide hi:त्वरित फुरिअर रूपान्तर id:Transformasi Fourier cepat it:Trasformata di Fourier veloce ja:高速フーリエ変換 ko:고속 푸리에 변환 nl:Fast Fourier transform pl:Szybka transformacja Fouriera pt:Transformada rápida de Fourier ru:Быстрое преобразование Фурье sr:Брза Фуријеова трансформација sv:Snabb fouriertransform ta:விரைவு ஃபூரியே உருமாற்றம் tr:Hızlı Fourier dönüşümü uk:Швидке перетворення Фур'є

Retrieved from "https://sound.bai-hua.org/index.php?title=快速傅立叶变换&oldid=122"