一步一步理解快速傅立葉變換

一. 傅立葉變換的原理

任何一種變換方法都是為了我們能夠更簡單的解決問題，傅立葉變換亦是如此。我們平時是如何去分解一個複雜的問題呢？一個經典的方法就是把這個複雜的問題分解成為多個簡單的可操作的子問題。

讓我詳細說來：

1.傅立葉級數

對於任何一個函數的曲線，傅立葉認為任何一個複雜的曲線都可以分解為正弦函數與餘弦函數的疊加！

函數的疊加： 

  

而對於傅立葉級數的思想我們可以把一些複雜的曲線分解成這樣： 

 
2.時域與頻域

剛才就說道傅立葉變換是一種變換，那到底由什麼變成了什麼呢？其實，傅立葉變換就是把時域中的函數轉換成了頻域中的函數，以方便我們做數據的處理，處理完之後，用同樣的原理再做一次傅立葉逆變換就得到了我們所需的結果！

域的意思是範圍也是一種視角，就如高中學過的定義域一樣。時域的意思就是隨著時間變化的視角觀察到的變化，可以理解為直角坐標系中x軸的意義為時間；頻域就是頻率域，意義也同樣如此。

對於上圖矩形波分解出來的函數們，他們在頻域中是這個樣子的： 

 

這兩者之間的轉換就是傅立葉變換！

3.再說傅立葉級數

還記得高中老師講過的三角函數麼？ 

 

動圖請看：File:Fourier series square wave

我們通常認為在時域中的基本單位是 1 ，而由上圖可以看到正弦波其實是圓在直線上的投影，所以我們可以認為在頻域中的基本單位是一個圓掃過的周期。

他們兩個的聯繫就是這樣的： 

 

從左面疊加看過去就是合成的時域圖像，而從右面看過去投射出的就是頻域圖像了。

這樣的好處就是用一個個的線段把每個波表示了出來，要是需要從複雜波中去掉某一個波的時候就十分方便了。

而在更深的層次，我們可以理解為把連續的函數抽象為了單個離散的點與線的集合！

4.連續傅立葉變換

在傅立葉級數中，我們把一個在時域連續的周期函數轉化成了在頻域離散的非周期函數，而連續傅立葉變換是把一個在時域連續的非周期函數轉化為一個頻域連續的非周期函數。（說著繞嘴，一個字一個字的讀） 

 

還記得傅立葉級數中分解的那幅圖麼： 

 

被分解的曲線密度越來越密，讓我想想一下，當曲線密集到一定程度的時候，一個曲線的左面無窮小處與右面無窮小處都存在另一個曲線的話，那麼我們認為他們之間是相互連續的，他們的疊加就構成了連續傅立葉變換，不過這時由於無窮小處也存在需疊加的曲線，那麼計算時的+號就應該變為積分符號了：

二. 離散傅立葉變換（DFT）

傅立葉變換存在4種變體，它們都是線性積分變換：傅立葉級數，連續傅立葉變換，離散傅立葉變換，離散時域傅立葉變換。

我們來看一下4種變換的圖樣與特性： 

大家都知道，計算機系統是離散的，所以離散傅立葉變換在計算機領域得到了很大的應用與改進。 
DFT的圖像上不是一條連續的曲線了，而是一個一個點組成的圖像，所以說它的公式是這樣的：

這裡也給出 逆離散傅立葉變換公式(IDFT)：

它最大的作用之一就是求卷積！

1.卷積

百度百科上講的卷積很抽象，這裡有一個「通俗「的例子：假如我打了你一拳，這一拳的痛感會在1個小時之內消失，但你這一個小時之內的感覺也不一樣。我們記使用最輕的力打你一拳時，你在這一個小時內的疼痛感覺圖像為2h（t），那麼用兩倍的力打你就是2h（t），3倍力3h（t）……f（n）倍力就為f（n）h（t）。如果我每秒打你一次，連續打你半分鐘，那麼你在一個小時零半分鐘之內都會感到疼（每一次有新痛感來臨時上一次痛感都會變化），那麼怎麼計算你在這一時間段內某一秒時的痛感呢？就是把之前的所有痛感疊加起來啊！Y(t) = 0 + …… + f（n）h（t-n），然而當你每次被打一拳的前無窮小時刻與後無窮小時刻都被打了一拳，那麼我就認為打你是一件連續的行為，這時候就要把加和改寫為積分了：

而卷積定理指出：函數卷積的傅立葉變換是函數傅立葉變換後的乘積！ 
也就是說要求兩個函數的卷積，我們就可以先把兩個函數做傅立葉變換，之後算出它們的乘積，再做一次傅立葉逆變換就得出結果了！大學計算瞬間變成小學算術啊有木有！

2.點值表示法

可以粗略的理解為卷積是兩個相關函數的乘積（這麼說並不準確）

求多項式乘法的本質就是求卷積，像這樣： 

可是如果單單就是這樣相乘就沒有什麼可優化的了，可是讓我們想一下多項式其實也是一個函數，每一個x值都對應一個y值，那麼就如兩個點確定一條直線一樣，我們也能用幾個點確定一個多項式，這就是點值表示法。 
如果多項式的最高次項的次數為n，那麼在常數已知的情況下，我們只需要知道n個不同點就能聯立求得多項式了，那麼我就可以用這n個點的集合表示這個多項式了。只要把點值對應相乘就可以了，複雜度為O（n）。 
那麼問題來了，如果單純的代值進去，運算量會非常大，複雜度為O（n²），這時候就需要用到快速傅立葉變換了。

三. 快速傅立葉變換（FFT）

快速傅立葉變換能在O（nlogn）時間內求出卷積的結果。

1.單位復根

FFT正是利用了單位復根作為旋轉因子把DFT二分才獲得了如此高的效率。

定義：若 ，則稱W為1的復根，也稱作單位復根。單位復根必有n個，它們是

為了方便們定義：

它的性質為以下這些：

Ok，基礎知識補充完成，看不太懂的話無所謂，接著往下看

2.FFT原理

回顧一下DFT的公式： 

當x[]數組存放的是複數，我們可以用如下代碼表示：

complex x[10],y[10]; for(int j = 0 ; j < n ; j++)        y[k] += x[j] * pow（e，-2πijk/n）; 
註：根據公式，調整一下
 
原理其實是一樣
首先讓我們把整個公式拆掉：所有數 = 奇數 + 偶數 （j = 2r + 2r+1）
再把把奇數項拆開：記x0[r] = x[2r] , x1[r] = x[2r+1]
計算量就從（n-1）變成了（n/2-1）瞬間減少了一半啊！
然後我們把它抽象出來就是：
仔細觀察上述公式，有沒有發現x0[r] 與x1[r]裡面的值為
所以：
這裡r = （0，…… ， n/2−1）
的確變快了但是這也就是說我們只能求出一半的值來，那後一半要怎麼計算？
回顧一下單位復根的折半引理：
我們就會發現當r=( n/2 , …… ，n-1 )時：
這兩半居然是對稱的！！
由於k對應著變換後的（0，……，n-1），我們就可以把FFT公式寫成：
這個就是快速傅立葉變換的原理了，這樣二分下去就可以在O（nlogn）複雜度下算出結果。
逆快速傅立葉變換同理，但要記得根據逆DFT，結果要除以n。
3.FFT的具體計算方法
FFT的計算方法被稱為蝶形計算：  
經過箭頭代表乘上箭頭的值，匯聚代表相加，這樣正好對應了我們的公式。
以8點FFT（8個數的）為例： 
它的運行順序有些奇怪是 0,4,2,6,1,5,3,7 ，看看第一個蝶形運算，我們把它們變形： 
0,2,4,6 
1,3,5,7 
兩個兩個看！是不是找到規律了，其實就是奇偶分開並且讓前一半的開頭與後一半的開頭對齊運算。 
我們把這種方法叫做倒序運算，這裡的倒敘指的是每個數按位倒敘： 
0 → 000 → 000 → 0 
1 → 001 → 100 → 4 
2 → 010 → 010 → 2 
3 → 011 → 110 → 6 
4 → 100 → 001 → 1 
5 → 101 → 101 → 5 
6 → 110 → 011 → 3 
7 → 111 → 111 → 7
但是要是把每個數按位逆序的話寫成程序計算量就會很大了，我們接著找規律： 
仔細看，是不是只有1,4與3,6互換了位置而其他都不變？ 
0與n-1這兩個數逆序與正序是相同的，這也限制了FFT中的n必須為2的n次方 。
對於中間的數（1~n-2），一組一組的看它們二進位的每一位，從低到高看，可以看出每一組是按高位值為0或者1 來決定分組的。我們設兩個指針i=1，i是可能被置換的數，一直把他歷遍到結尾，它的二進位是000……1；j = n/2 ，j是可能要置換的值，我們通過一系列的操作去維護它，它的二進位恰好是1……000（不信試試看）。
接著往下找，以後的j 要比i大才有意義，剛才看到了一個數加上n/2也就是可以把最高位由0→1 ， 同理加n/4是次高位由0→1，加n/8是次次高位由0→1……，so，我們每一次從右往左尋找，遇1變0，遇0跳出並且把該位變成1，這樣當i<j 的時候置換他們，否則就接著找，這種方法叫做叫做二進位原地逆序置換。 
讓我們試一下16點的離散傅立葉變換： 
00 → 0000 → 0000 → 0 
01 → 0001 → 1000 → 8 
02 → 0010 → 0100 → 4 
03 → 0011 → 1100 → 12 
04 → 0100 → 0010 → 2 
05 → 0101 → 1010 → 10 
06 → 0110 → 0110 → 6 
07 → 0111 → 1110 → 14 
08 → 1000 → 0001 → 1 
09 → 1001 → 1001 → 9 
10 → 1010 → 0101 → 5 
11 → 1011 → 1101 → 13 
12 → 1100 → 0011 → 3 
13 → 1101 → 1011 → 11 
14 → 1110 → 0111 → 7 
15 → 1111 → 1111 → 15
其實是一樣的！
這段程序可以這樣寫： 
先給出複數結構體
進行倒序運算：
倒序之後我們開始按照蝶形計算的方法計算：
這樣我們就能快速的對函數進行離散傅立葉變換了！
4.快速傅立葉變換到底有多快？
粗略認為最高次數相同時： 
DFT計算時多項式對應相乘，複雜度為：變換O（n²）+ O（n²）,頻域相乘O(n)，逆運算O（n²），總複雜度為O（n²）； 
FFT計算時，變換O（nlogn）+ O（nlogn）,頻域相乘O(n)，逆運算O（nlogn），總複雜度為O（nlogn）； 
作出圖來是這樣： 
 

當1024點時FFT的計算量大約是DFT的1%！
四. 複習一下
傅立葉變換把函數從時域變換到頻域
時域卷積對應著頻域乘積
點值表示法把函數離散化了
離散傅立葉變換可以把點值變換到頻域
點值在頻域的乘積結果，逆離散傅立葉變換後得到時域函數的卷積結果
快速傅立葉變換巧妙利用了計算順序使計算結果重複利用減少了計算次數
五. 在ACM中的應用
超大整數相乘
大數運算講解題集中的題目
2.組合數學計數
直接看例子講：hdu4609 
題意：給出n個邊的長度，求任取3條能組成三角形的概率
題中的例子a[4] = {1,3,3,4} 把每個長度的邊的個數記下來 num[5]={0,1,0,2,1}， 
長度為3的邊有兩條，為2的有0條。然後把它們兩兩組合起來{0,1,0,2,1}×{0,1,0,2,1}， 
看出來了麼，是不是就是求卷積！得出的結果num[4*2+1]={0,0,1,0,4,2,4,4,1},它的意思也是計數個數組合後長度為0 的有0 個，長度為1的有0個，……，長度為7的有4個，為8的有1個。
為了計算方便，我們把輸入的數組a排序，把每一個邊都當作最長邊去算，這樣的話雖然簡便了但還是會有重複的值，把他們刪去就得到最終的結果了。 
具體的解釋邊看代碼邊說（已AC）：
本文轉載自csdn博客，版權歸原作者所有，代碼未測，不知是否有效。