tucao

先要吐槽一下簡書，居然用不了MathJax囊咏，導(dǎo)致所有的公式要寫成Latex亿笤，真是反人類
||Φ|(|T|Д|T|)|Φ|||

void 吐槽()
{
最近剛看完FFT，頓時心中真是成百上千的mmp符匾。倒不是方法本身的問題，而是現(xiàn)在的書籍和資料很多真是對萌新極其不友好瘩例，尤其是國內(nèi)一些所謂的“入門到入坑”系列啊胶，總是把一些簡單的問題說復(fù)雜，而且描述同一個公式還版本多種多樣垛贤，并且把一些應(yīng)該放到后面的知識排到了前面焰坪，要么讓你望而生畏，要么間接up你查閱資料的能力聘惦。
好的某饰，吐槽完畢。
return
}

這篇簡書只是關(guān)于DFT筆記的匯總，大部分是個人主觀的理解黔漂，對DFT理解不深诫尽，文中一些不嚴(yán)謹(jǐn)?shù)牡胤綒g迎指出糾正。

開篇之前要先感謝一位師兄炬守，將BASIC版本的Bit Reversal Sorting 翻譯成了C版本牧嫉，本文FFT的程序才能完整實現(xiàn)。
你部長還是你部長 <)减途。(>

目錄

• DFT基本理論
• DFT四種算法
• 使用DFT去噪聲與解調(diào)

DFT基本理論

DFT全名Discrete Fourier Transform (離散傅里葉變換) 酣藻，是傅氏變換家族中用于處理在橫縱軸上數(shù)值離散且具有周期性的信號。相比于另外三個鳍置，這種變換方法處理有限并且周期延拓的信號辽剧，所以適用于計算機處理。

FT主要的目的是把時域信號中的各頻率成分的振幅和相位分別提取出來税产，從而將時域信號轉(zhuǎn)換位頻域信號中的幅頻譜和相位譜上怕轿。
對于一維的時域信號，要得到頻譜只需進行一次DFT砖第；而對于二維（圖片）以及二維以上的信號撤卢，需要進行多次的不同方向上的DFT，本文只涉及一維信號的DFT梧兼，對于變換圖片的頻譜請參考相關(guān)資料并類推即可。

DFT四種算法

對時域信號進行DFT的方法有很多智听，簡wo單zhi的hui有以下幾種：

一羽杰、聯(lián)立方程組求解 Σ( ° △ °|||)︴

根據(jù)傅里葉分解的原理，任何信號都可以多組正交的正弦曲線疊加而成
每組正弦曲線的頻率ω不同到推，相位Ф不同考赛，振幅A不同
根據(jù)香農(nóng)采樣定理，采樣率應(yīng)為信號最高頻率的兩倍以上莉测，取最低要求情況颜骤，若剛好位兩倍的情況，ω的范圍是0（直流）到0.5f（采樣頻率的一半）捣卤。
列出方程組求解得到多組的解A和Ф忍抽，求得幅頻數(shù)組Mag[]和相位數(shù)組Phase[]。

線性方程組

數(shù)學(xué)專業(yè)的同學(xué)了解一下……
嗯董朝，1024個點為1024個線性方程組鸠项，解吧:-)

二、相關(guān)性求解法（傳統(tǒng)方法）

方法一是從數(shù)學(xué)的角度來理解DFT子姜，可以作為理論上的用法祟绊，但是由于計算量過大，不能用于計算機處理。

若從求相關(guān)性的方法來計算牧抽，則可用于計算機處理DFT嘉熊。
相關(guān)性求解的方法屬于一種類卷積的運算，通過將原始信號與目標(biāo)頻率信號求相關(guān)扬舒，可以得到二者的相關(guān)度记舆，再求代數(shù)和，即可得到該頻率下的信號相關(guān)度呼巴。

求相關(guān)

這里改方法根據(jù)原始信號的種類不同又可以分為兩種：

• 實信號處理模式（更快泽腮，但是公式混亂）
  對于N點的實數(shù)信號，根據(jù)以上的方法衣赶，可以分解為N/2點的頻域數(shù)組Re[]和Im[]诊赊。所以求相關(guān)性時應(yīng)該使用對應(yīng)頻率的正弦函數(shù)sin和相同頻率的余弦函數(shù)cos，求得的相關(guān)值分別放入Im[]和Re[]府瞄。

• 復(fù)信號處理模式（慢一些碧磅，但是統(tǒng)一度好）
  對于N點的信號，將其當(dāng)作復(fù)數(shù)處理遵馆，相當(dāng)于有N點的實部和N點的虛部（對于純實數(shù)信號鲸郊，N點虛部全為0），再將其與旋轉(zhuǎn)因子（有的版本稱相位因子）W_N

進行相關(guān)度運算货邓，得到實部Re和虛部Im秆撮。
值得注意的是，此時Re數(shù)組和Im數(shù)組長度均為N而不是N/2换况。

DFT&IDFT

（0 ~ N/2）部分的數(shù)值和“實數(shù)模式中得到的結(jié)果一摸一樣职辨，稱之為正頻率，而（N/2 ~ N）的數(shù)值為負(fù)頻率戈二。
負(fù)頻率由于本身沒有物理意義舒裤，且和正頻率水平鏡像對稱，為冗余信息觉吭，故可以舍棄不運算加快速度（但是負(fù)頻率在FFT中有著其存在意義）腾供。

得到相關(guān)性數(shù)組Re[]和Im[]，是直角坐標(biāo)系下的頻譜鲜滩，在運算上有作用伴鳖，但是不直觀。

極坐標(biāo)形式的幅頻圖

根據(jù)高中學(xué)過的Asin(ωt+Ф)=Xsin（ωt）+Ycos（ωt）可知绒北，Re[]和Im[]的信息可以等價保存為幅頻數(shù)組Mag[]和相位數(shù)組Phase[]黎侈，將直角坐標(biāo)系轉(zhuǎn)化為極坐標(biāo)系，使得數(shù)據(jù)更加直觀（尤其是頻譜部分）闷游。

三峻汉、對稱性求解法

對稱性方法比較容易理解贴汪，但是計算量卻十分的大（僅次于方程組求解法）

我們知道，頻域的一個δ(n)沖激響應(yīng)（一個點）相當(dāng)于時域的一個響應(yīng)頻率的正弦曲線休吠，根據(jù)時域與頻域的對稱性扳埂，時域的一個δ(n)也相當(dāng)于頻域的一條正弦曲線
（算是是自然界中最神奇的對稱之一）
所以，可以從時域逐點求得對應(yīng)的正弦曲線瘤礁，在頻域進行累加運算阳懂。
時域上每個δ沖擊函數(shù)都有位移，相當(dāng)于δ(n)與Δ(n+1)進行了一次卷積柜思，所以頻域的每次累加要先乘于一個對應(yīng)頻率的正弦曲線岩调。（因為Δ(n)的傅里葉變換對是sin）。
說白了就是左邊一個個數(shù)幅度赡盘，然后右邊慢慢加

四号枕、基2法（“饑餓法”！陨享？）

還有一種叫基4法（這篇比較全面）：
http://www.doc88.com/p-7979920042436.html

基2法包括DIT時域抽取和DIF頻域抽取葱淳，這里只簡單介紹以下DIT。
基2法全名是 DIT-FFT-基2 算法抛姑。
此處FFT是Fast Fourier Transform即快速傅里葉算法赞厕，F(xiàn)FT有很多種，需要注意的是FFT不是新的傅里葉分解方法定硝，它只是DFT的一種方法改進皿桑，就是公式還是原來的，只是優(yōu)化了喷斋。
相比于傳統(tǒng)的DFT唁毒，F(xiàn)FT大約能快個300倍左右，并且數(shù)據(jù)量越大星爪，F(xiàn)FT越有優(yōu)勢。

• 傳統(tǒng)的DFT為：

DFT

對于N點的時域信號粉私，算一次DFT顽腾，需要進行 _N2 次復(fù)數(shù)乘法和累加運算，這就是傳統(tǒng)DFT的局限性之一诺核。
在100Mz的處理器上抄肖，計算 ₂10 個采樣點，即便預(yù)先算好旋轉(zhuǎn)因子制成一個數(shù)組表窖杀，估計也要十幾秒漓摩，那么對于更大的點的速度將會是不可接受的（Time out啦兄dei）。

之所以這么慢入客，是因為傳統(tǒng)DFT的計算有著大量的冗余管毙。利用旋轉(zhuǎn)因子W_N^nk本身的對稱性腿椎、周期性、可約性以及特殊值夭咬，簡化DFT的計算步驟得到FFT算法啃炸。

• 分而治之
  FFT算法利用的是“分而治之”的算法思想，對于DIT基2型的FFT卓舵，將原始時域信號X[k]一分為二南用，一邊是偶次點X₀[k]，另一邊是奇次點X₁[k]掏湾。

則DFT(X[K])=DFT(X₀[k])+DFT(X₁[k])

根據(jù)旋轉(zhuǎn)因子的特性裹虫，將公式變形后可以得到

X[k]=F₀[k]+(W_N^nk)F₁[k] k=0,1,2,……,N-1

二點運算

而F[k]又可以同樣地往下拆，直到只剩下一個點融击。
然而一個點的時域數(shù)值就是頻域數(shù)值筑公，由此回溯，最終得到Re[]和Im[]砚嘴。

蝶形算法

使用DFT去噪聲與解調(diào)

DFT的應(yīng)用很多
不如說十酣，DSP就是圍繞著傅里葉變換而發(fā)展起來的

例如，語音通話去背景噪音际长、音樂人聲去除耸采、語音識別的聲學(xué)特征提取、雷達索敵工育、聲納探測虾宇、神經(jīng)網(wǎng)絡(luò)……

多倒是挺多的，不過要做成實用的還有很多邊沿技術(shù)和交叉學(xué)科要學(xué)

看了FFT之后興沖沖地想去親手寫一次變換程序如绸，結(jié)果還是太年輕了嘱朽，卡在了第一步的位反轉(zhuǎn)分類算法，寫了幾個版本就是得不到正確的結(jié)果怔接。
BASIC版本的又太混亂看不下去搪泳，無奈向部長求助，結(jié)果是一直把1(0001)當(dāng)成F(1111)了......
(;′??Д??`)

這次就不用Matlab（MatLab功能太強扼脐，一些細(xì)節(jié)問題都幫我們處理了岸军，就不知道問題出在哪了），然后 .NET 又不會瓦侮，Qt又不會艰赞，GTK+也不會
（那我到底會什么(ノへ￣、)肚吏。方妖。。罚攀。党觅。雌澄。）
好吧，會一點Java仔役，那就用swing吧

code

利用一個int數(shù)組（double也行掷伙，不過Panel的繪圖是像素位單位，還要再轉(zhuǎn)回去比較麻煩）又兵，保存1024個抽樣點的時域值任柜，再創(chuàng)建ReX[]和ImX[]保存直角坐標(biāo)系下的頻域值，MagX[]和PhaseX[]保存極坐標(biāo)系下的頻域值沛厨。

需要注意的幾個點（坑）：

• IDFT逆離散傅里葉變換之后要乘一個-1宙地，不然圖是反的，至今不知道為什么
• 可以定義一個PI常量為3.1415926
• 若ReX[k]為0逆皮，則不能求arctan（ReX[k]/ImX[k]）宅粥，要根據(jù)ImX[k]的正負(fù)判斷相位
• 當(dāng)MagX[k]小到可以忽略時，相位Phase[k]無意義

image.png

窗口左邊的按鈕是產(chǎn)生時域波形电谣，使用sin函數(shù)生成一個1024長度的數(shù)組秽梅，右邊按鈕是FFT，將波形轉(zhuǎn)化為頻域信號剿牺。

--------------------------分割線----------------------------

先來一個sin波形企垦，頻率為（1/512）Hz

F=(1/512)

FFT后：

FFT

雖然不是很明顯，但是最左邊有一個高高立起的（flag）頻率為一的δ沖激信號晒来，其他地方都是0钞诡，說明波形中只有一種頻率。
（右邊一坨的是相位湃崩，用于顯示突變的點）

換成f=(10/512)后：

F=(10/512)

FFT

頻譜往右挪動了一丟丟荧降，但還不是很明顯

F=350/512：

F=(350/512)

FFT

很好，已經(jīng)確定DFT分析成功攒读，但是現(xiàn)在有什么用呢朵诫？
對于一種波形DFT確實沒什么好做的，但是現(xiàn)實中的波形肯定不止一種薄扁，而是多種波形混合而成

比如：
f=3/512混合上f=50/512混合上f=150/512的情況是：

混合信號

假如f=3/512是語音信號（當(dāng)然不可能這么低泌辫，只是舉例說明）
f=50/512是載波信號（語音信號的載體）
f=150/512是噪音信號

現(xiàn)在三者一起傳過去，怎么樣從上面混亂的信號里得到我們想要的聲音呢九默？
不妨先用傅里葉分析看一下頻譜：

FFT

for(int i=140;i<MagX.length;i++)
{
MagX[i]=0;
}

純信號

現(xiàn)在是純凈信號墨礁，比之前好看很多

要提取我們的聲音信號只需把f=50/512的載波去除就好了幢竹，方法和上面一樣：

for(int i=40;i<MagX.length;i++)
{
MagX[i]=0;
}

語音信號

顯然有一些失真，這是因為在FFT的時候沒有給原時域信號加窗函數(shù)平滑處理恩静，導(dǎo)致時域分幀后能量泄露焕毫，引起邊緣混疊，再變回來就變成這樣驶乾。
可以選擇sinc窗函數(shù)或者漢明窗函數(shù)：
再變回來就好一些了：

平滑處理

總結(jié)

程序員的優(yōu)點除了懶惰外邑飒，最大的就是動手的機會多，DFT看似容易级乐，其實很多問題要在實際編程中才能發(fā)現(xiàn)疙咸，解決細(xì)節(jié)問題。

傅里葉變換是DSP的核心风科，一切從這里出發(fā)撒轮，一切又回到這里。

之前在某乎上看到一句話寫的很好：

時域上的變化無常的錯綜復(fù)雜的曲線贼穆，只是頻域上的一個點
你眼中看似落葉紛飛變化無常的世界题山，實際只是躺在上帝懷中一份早已譜好的樂章。

最后上一張傅里葉的微笑：

B in your heart