因為要移植CSK得寫快速傅立葉變換的演算法，還是二維的，以前在pc平臺上只需呼叫庫就可以了，只是有點印象原訊號和變換之後代表的是什麼，但是對於離散傅立葉變換的來龍去脈忘得已經差不多了，最近要用到，於是重新來學習一遍，翻出了自己大三當時錄的吳鎮揚老師講的數字訊號處理的視訊，DFT-FFT這裡老師講了有10講之多，但每講都不是很長，20分鐘左右，這裡記錄一下學習的過程，前面的推導有點多，簡書又打不了公式，mathtype的直接複製也不過來，截圖又太麻煩，也為了自己再推導一遍，手寫了前面一部分的內容。圖片形式傳上來。
簡單說幾句：DTFT有了之後為什麼還要搞出來一個DFT呢，其根本原因就是因為DTFT的頻域是連續的，無法用計算機進行處理。根據我們之前得到的的傅立葉變換的規律：

連續----非週期
離散----週期
週期----離散
非週期----連續

根據這四個規律，只有一種情況可以保證兩個域都是離散的，那麼就是週期離散的訊號，對應的另外一個域是離散週期的，但是我們又知道週期的訊號是沒有傅立葉變換的，因為在整個Z平面是不收斂的，這就是一個問題了，討論就從這裡開始。

從DTFT到DFS

從DFS到DFT

簡單的來說，DFT是針對有限長序列的，那麼怎麼來做DFT呢，這裡的做法是找到其對應的週期延拓序列，做DFS，然後再擷取主值序列。相當於是把DFS引申到有限長序列裡來，之所以能夠引申是因為DFS的時域和頻域都是週期的，而且週期還是相同的。

上面討論了DFT的三個性質，分別是線性，迴圈移位和迴圈卷積,關於迴圈移位和迴圈卷積有必要說幾句：

DFT的迴圈移位和迴圈卷積分別對應著DFS的線性移位和週期卷積，這兩者實際上是有著很強烈的關係的，甚至在某種情況下是完全一樣的：那就是當我們只關注DFS的一個週期時，迴圈卷積和線性卷積是一樣的。(所謂迴圈移位就是從一端移出去的要從另一端移進來)
這裡的根本原因在於：當我們只關注一個週期時，週期序列的線性移位和和非週期序列的迴圈移位的結果是完全相同的。

看個例子：

移位

上面都是向右移動兩個單位，如果只關注主值的話，迴圈移位和線性移位的結果是完全一樣的。正是由於這個原因，所以才有了迴圈卷積，卷積裡也有移位的操作，迴圈卷積和週期卷積來說，如果只關注一個主值序列的話，那麼結果是一樣的。
關於迴圈卷積還有一些更重要的性質。

迴圈卷積（針對兩個有限長序列）

1. 若：
   

   頻域點積

   則有：

   
   
   

   時域迴圈卷積
   
   主要意思：時域迴圈卷積相當於頻域點積。
2. 計算過程：
   ① 有限長序列構造週期序列
   ② 計算週期卷積
   ③ 週期卷積取主值
   迴圈卷積是可以用DFT來做的，因為1的性質，如果直接用DFT來做的話，計算量實際上是要比卷積大的， 但是因為我們有FFT，所以許多情況下迴圈卷積還是可以通過FFT來做的，而且計算量會大大減小，特別是其中的一個序列不變的情況下（比如濾波器，那麼這個濾波器只需要做一次FFT）。
   這樣可能又有疑問了：對於線性系統來說（比如脈衝響應），是通過線性卷積來算的，移出去就不會再移回來了，那麼貌似迴圈卷積又沒什麼用了？答案肯定是否定的！，如果是這樣我們討論著半天不是白說了？實際上，在某些情況下是可以用迴圈卷積來計算線性卷積的，下面討論這種情況，看下要滿足什麼條件。
3. 有限長序列的迴圈卷積和線性卷積。
   上面說了：迴圈卷積在一定條件下可以計算線性卷積。
   假設有兩個有限長序列：x(n)和y(n)，長度分別為N,M，那麼線性卷積的長度是N+M-1（算一下就知道了）。
   那麼我們如果對這兩個序列做迴圈卷積呢？要做迴圈卷積，序列長度首先得一樣，那麼怎麼變得一樣呢？在後面添0。添多少呢？現在還不知道。
   我們假設添到L(L>max(M,N))，添0對線性卷積是沒有影響的，為了分析兩者的迴圈卷積，先看其週期延拓：
   
   

   週期延拓
   
   那我們先把週期卷積算一下：
   
   

   週期卷積
   
   結論就是：週期卷積相當於線性卷積(上面紅色的字打錯了，我不小心把mathtype關了，懶得改了)的週期搬移。而迴圈卷積是週期卷積的主值序列。
   這樣就說明說明呢？說明週期卷積是線性卷積的週期延拓的關係，延拓的週期是L。
   這時候就要關注混疊了，因為L必須足夠長才能保證搬移的時候不會發生混疊，結合上面線性卷積的長度，那麼L的長度最少就是L>=N+M-1，這樣才不會產生混疊。這樣取主值區間才能取到線性卷積的結果。
   這樣就可以用迴圈卷積做線性卷積。

共軛對稱性

設x(n)的共軛復序列是x*(n).則：
DFT[x*(n)]=X*(N-K)
證明：

由於W是週期的，且週期是N，所以可以寫作：
看這個結果和DFS其實是一樣的，這裡只不過把它移動到主值區間上罷了。
分別x(n)看實部和虛部：

根據線性變換的性質，實部和虛部分別做DFT，那麼可以得到：

同理：

那麼顯然是有：
為什麼把上面的叫做奇對稱和偶對稱分量呢？看下面的推導，對於偶對稱來說：
即：

同理的：共軛奇對稱分量為：

實際上：共軛偶對稱就是實部DFT的結果，共軛奇對稱是虛部DFT的結果，滿足這樣的關係。
那麼對於純實數序列來說，其只存在共軛偶對稱部分，表明實數序列的DFT滿足共軛偶對稱性，利用這一特性，只要知道一半數目的X(k)，就可以得到另一半，這一特點可以在DFT運算中加以利用，提高運算效率。。

還有一個例子：對於一個實數序列（長度為n）來說，其DFT是2n個點（那個虛數），資料量是增加了一倍了的，而他們之間只是線性變換而已，為什麼多了這麼多資料？實際上用共軛偶對稱可以很好解釋，實際上沒有多，有一半的資料就可以通過共軛對稱性得到另外一半。

另外值得一提的是：DFT的兩端是對稱的，即x(n)與X(K)是對稱的，那麼如果把頻域作為原始資料，依然可以找到時域的共軛偶對稱和奇對稱部分。不細說了。

選頻性

DFT具有選頻性。

對於復指數函式：
進行取樣：

其中
其中q是整數，w0=2pi/N時：

其傅立葉變換為：

這個式子前面討論過，這是n個模為1的向量，只有k=q+sN時是1，其他時候都是0，在這個主值序列裡就是k=q。

這個是顯而易見的，如果輸入序列只有一個頻率，那麼對這個序列取樣再進行DFT就應該只有一個頻率是有值的。
當如數頻率是qw0時，變換X(k)的N個值中只有X(q)=N,其餘均是0，如果輸入訊號為若干不同頻率的訊號的組合，經離散傅立葉變換後，這些頻率對應的X(k)應有對應輸出，因此，離散傅立葉變換演算法實質上對頻率有選擇性。

看個題：

一個餘弦序列，頻率是qw0,這樣得到的頻譜就是兩條線，一個q，一個其對稱位置。這也驗證了剛才共軛偶對稱的性質。

DFT和z變換。

DFT是DTFT頻域取樣的結果。

這樣取樣也會有問題的，就相當於通過一個柵欄來觀察DTFT的結果，肯定有資訊丟失，可能就會有有問題，到底取樣密度怎樣才是好的。
實際上從上面討論中可以看到了，N個點的DFT是N個點，所以在頻域的取樣數M的條件應該是：M>=N。
證明：
這個證明的關鍵也是在於交換了一個積分順序，結論是很直觀的，就是在頻域取樣就相當於是在時域週期延拓，那麼對於一個長度為N的序列，只有當以大於N的週期去延拓時才不會發生混疊，所以頻域的取樣應該至少是N個點。把這個直觀的結果記住，具體的證明如果要用到查書即可。

中間吳老師還講了帕斯瓦爾定力，一些訊號處理的流程，把整個知識串起來的，我快快得看了一遍。已經等不及去看FFT了。
頻譜洩露這一次才真正理解了，頻譜洩露就是加窗時發生的，離散週期訊號要進行DFT時要進行截斷，如果不是整週期截斷，做DFT得到的頻譜就會發生洩露，本質的原因就是週期延拓的時候就不是原先的訊號了（因為沒有整週期截斷），這個是個充要條件。
其他的就不說了。

從DFT到FFT

DFT並不是新的演算法，但是直到FFT的發現，才讓DFT真正運用到工業和生活中，1965年cooley(IBM)和Tukey(MIT)提出了2FFT（2的冪次）演算法。這樣使得DFT的計算量提高了1,2數量級（與N有關）。這個是基2的DFT演算法。
實際上還有其他的演算法，一個典型的演算法是：76年IBM的S.Winograd博士推導了WFTA演算法，使FFT得運算再次下降很多，但是由於其用到了取模運算來對映地址，而這種定址是非常耗時的，所有大資料的FFT都不用了，只在一些小範圍內還有運用（查表方式取模）,這段是吳老師科普的。
下面主要介紹基2的FFT的演算法：

DFT的計算。

首先我們看下要進行n點DFT運算時要進行的計算量：

實際上這兩者變換隻是差了一個指數的負號和一個常數，其計算量是完全相同的。這裡只討論正變換的計算量：
一般來說，我們認為x(n),W,X(k)都是複數，那麼每計算一個k，需要N的複數乘法，以及N-1個複數加法（N個相加），一共需要計算N各k，所以一次N點的DFT要用到：
N*N的複數乘法和N(N-1)的複數加法。
而實際上覆數的乘法和加法都是通過實數完成的，每一次複數乘法需要4次實數乘法和2次實數加法，所以最後換算到實數計算量：

每一次N點的DFT所需的計算量：4N^2次乘法和2N(2N-1)次加法。
DFT有幾個特點：

1. 係數W是週期的和對稱的。

   
   
   
2. 計算量是正比於N^2的，所以當N比較小的時候，計算量是可以接受的，所以一個思路就是把大的序列拆成小的序列。

基2的FFT主要有兩種，按時間抽取的和按頻率抽取的。重點說按時間抽取的，即按照n抽取的。

按時間抽取的基2FFT

首先我們假定序列x(n)長度是2^M,不夠的話補零。

首先將序列分成奇偶兩組：
帶入求DFT得公式之中：

由於：
則有：
這兩個求和實際上還是DFT，分別是偶數序列和奇數序列的，點數均是N/2。然後加權求和。

其實我一開始很糾結這塊關於括號裡的2r，這個其實不要被表面矇騙了，雖然是2r,但是在這個序列裡還是代表的是第r個數，所有求和符號與W裡都化簡成了r，都是從0開始到N/2的自然數。

根據上面分析的，則有：

其中兩個x分別是偶數序列和技術序列，這樣就很明顯了，就是兩個分離的DFT。
但是我們知道，N點對應的DFT的點數是N點，這兩個都是N/2的，那麼對應變換再求和還是隻有N/2個點啊，但是X(k)應該是有N個點的，怎麼會少了一半呢？帶著這個疑問再往下看吧。

我們把上面的式子再簡寫一下：
其中G(k)和H(k)分別代表偶數序列和奇數序列的DFT。都是以N/2為週期的。
所以：

前面證明過：
所以有：
所以N點的DFT可以拆成奇偶兩部分，然後分別以不同形式加權（加和減）獲得，前一半的k做加法，後一半做減法，這樣就能得到N點的DFT結果。
如果我們簡化用a，b，w分別表示上面的式子，那麼我們需要計算兩個式子：

畫的這個圖就是簡化的地方：右邊的計算只需要一次複數乘法（wb只需要計算一次）和兩次複數加法。這就是一個蝶形運算。
到這裡歸納一下：一個偶數點的DFT可以分成奇偶兩部分，然後通過蝶形運算加權起來，我們這裡取得N是2的m次冪，所以這樣我們可以按照這個思想一致分下去，一直分到2為止。
先看一個簡單的：
以N=8為例，做一次奇偶分解，然後通過蝶形運算組合起來：

這個很好理解吧，還可以再拆一次，拆的過程省略了，直接看計算的結果：

x的下標是分解的結果，x1是第一次分解的偶數，x2是第一次分解的奇數，x3是x1分解的偶數，x4是x1分解的奇數，以此類推。
稍微有疑問的一點可能是做完N/4的DFT之後的因子為什麼是W(N-0)和W(N-2）,這是因為：
這樣就很清楚了。這樣表示是把所有的W因子都用N為底的來表示。寫起來好看也好算
2點的DFT本身就是一個蝶形運算了，相當於再分解了一次：

這就是N等於8時按時間抽取FFT的運算流程圖。
所以2的整數次方的DFT完全可以由蝶形運算計算，這樣就大大降低了計算量。
現在就只剩下一個問題了：如何分奇偶？如果要通過除以2判斷餘數的話，那麼資料量大的時候這種演算法效率還是不高的。
下面介紹時間抽取FFT的四個特點，根據這些特點才可以設計出高效的計算機計算演算法。

1. 蝶形運算
   

   2^m點的資料，共有M排蝶形運算，有N/2個蝶形（每兩個資料一次蝶形），每一次蝶形1次複數乘法，2次加法，m是N以2為底的對數，所以N點的計算量為：
   
   
   這裡寫的都是複數的運算量。
2. 原位計算
   這是很好理解的，這是個單向的運算，一旦算到第二排，第一排所有的資料就不再需要了，這樣的好處是一個流水操作，跟前面的是沒有關係的，可以節省了很多的運算空間。

3. 序數重排（亂序）
   在進行蝶形運算之前，是要把原始自然順序的資料要重新排列，這個如果沒有規律耗時的話這個演算法就不會成功了。
   值得慶幸的是，基2fft的這種亂序是有規律的。以這八點為例，我們把這種亂序的序號用三位二進位制表示，然後再中心對稱，看看結果是什麼：

   
   
   
   
   

   看最右邊的一列是什麼？就是序列的順序標號，這種操作顯然是可逆的。也就是說我把自然順序的地址反過來就可以得到新的地址，這種地址的規律性可以稱作按位翻轉的地址。
   一些專門為訊號處理的處理器內建了按位翻轉單元，能夠以更高的效率來進行fft運算。
   實際上這個規律不是說硬是這樣發現的，而是隱含在蝶形運算裡的。第一次分奇偶的時候是按照最低位分的，分好的奇偶在各自的序列裡同樣是按照奇偶來分的。這樣就是一位一位來分，反過來就相當於是從高位來。

4. 蝶形的型別成倍增加
   蝶形的型別主要區別在W上，取數的間隔也是這樣的一個規律。

這樣的話，有了這四個規律，基2FFT的整個計算就很明瞭了。不夠2的基數是可以添0處理的，這樣還能改進柵欄效應。

關於基於頻率抽取的演算法沒有仔細關注，只是簡單的看了一下，感覺和按照時間抽取整個過程是反過來的。計算量也是完全一樣的，時間順序不用重排，但是頻譜的順序是亂序的，依然需要重排。而且亂序的規律都是一樣的。原位計算也是滿足的，蝶形的型別是隨著次數的增加而減少。

還有一個更有意思的，如果把這整個計算系統當做一個系統的話，把輸入當輸出，輸出當輸入，當中的箭頭全部反向，這就是同一個東西了。我找了一張圖，可以看下。

image.png

吳老師講：這個在訊號流圖中被稱作轉置定理，如果有興趣的話可以去看看這個，我這裡不細究了。

當然還是有其他的演算法，基4的，N是組合數的，如果有興趣也可以找來研究，我瞭解到這裡就足夠了。

---

總結：至此為止，從DTFT開始，如何一步一步得來到DFT以及怎樣得到FFT的演算法，我覺得已經總結得很清楚了，中間有大量的公式都是在mathtype上敲好然後截圖過來的。還是花費了一定的心力，這種經典的理論如果看進去了就真的和一本好書一樣，常看常新，經常都會有新的理解。下面應該還會寫程式設計實現和二維的怎麼來算。如果寫了連結貼在這裡。

---

從一維到二維

本來想重寫一篇的，後來發現從一維到二維的推導是如此的明瞭和簡單，就放在這裡了：
訊號中的fft大都是一維的，影象是二維訊號，在影象中的頻譜分析都是一維的，所以有必要對二維的DFT進行分析。在看這個之前，應該要對一維的DFT有充分的瞭解，可以看這裡。
二維DFT的公式是怎麼來的就不仔細推導了，只是看著公式推導一下和一維DFT之間的關係。

二維DFT

對於M行N列的矩陣，我們定義其DFT為：

我們按照其求和順序寫開：
當依照求和順序把不相關的提出去的時候就很容易發現，這實際上是可以分解成兩組DFT的，可以對每一列先進行DFT，然後得到的結果還是這麼大的矩陣，然後在在此基礎上進行行DFT，這樣得到的最終結果就是二維矩陣的DFT，二維DFT我們也是依照這個思路去算的，DSP的函式庫裡提供了一維的DFT運算函式，應該是效率比較高的，可以去借助這個實現二維離散傅立葉變換。

---

1-18，修改了些錯別字。