摘自21IC
只要是理工科畢業(yè)的朋友，都學過傅立葉級數與傅立葉變換，但真正要與實際應用聯(lián)系起來，用它來闡述應用中的各類問題，我們總會感覺概念模糊，似懂非懂，不知從何說起。是的，作者和你一樣，常常有這樣的體會。現在，讓我與你一起重新學習傅立葉的基本理論和應用，最后還給出一份FFT（快速傅立葉變換）的源碼（基于C）。希望對你有所幫助。Let’s go！
1． 歷史回顧
談傅立葉變換，不能不說三角函數。三角函數起源于18世紀，主要是與簡諧振動的研究有關。當時的科學家傅立葉對三角函數作了深入研究，并用三角級數解決了很多熱傳導的問題。三角函數的展開式如下：
f(t) = (1/2a0) + (a1·cos(x)+b1·sin(x)) + (a2·cos(2x)+b2·sin(2x)) + …
其中，系數a和b表示不同頻率階數下的幅度。
成立條件：
n 周期性條件，也就是說f(x)描述的波形必須每隔一段時間周期T就會重復出現；
n Dirichlet條件，周期T內，有限的最大最小值，有限的不連續(xù)點；
任何區(qū)間內絕對可積；
研究目的：
把一個基于時間變量t的函數展開成傅立葉級數的目的是分解為不同的頻率分量,以便進行各種濾波算法。這些基本的組成部分是正弦函數SIN（nt）和余弦函數COS(nt)。
應用領域：
l 信號分析，包括濾波、數據壓縮、電力系統(tǒng)的監(jiān)控等；
l 研究偏微分方程，比如求解熱力學方程的解時，把f(t)展開為三角級數最為關鍵。
l 概率與統(tǒng)計，量子力學等學科。
2． 傅立葉變換
H(w) = ∫h(t)·e^jwt·dt, （區(qū)間：-∽～+∽，w = 2πf）
討論：這里為什么會選擇復指數的形式而沒有用正弦余弦表示？
答案：歐拉公式的引入使得這條經典的數學公式變得更簡單，即e^jx = cos(x) + jsin(x)
3． 快速傅立葉變換（FFT）
常規(guī)的傅立葉變換算法并不適用于嵌入式控制系統(tǒng)，原因是運算量太大（涉及到復數運算），比如離散的傅立葉變換等同于用序列Y(n×1列矢量)乘以n×n矩陣Fn，需要n×n次乘法。若n=1024,則是104，8576次乘法運算。哇，這么多呀！什么概念呢？如果你選用的CPU單周期指令為25ns, 單周期也可以完成一次乘法運算，那么要計算1024點的傅立葉變換則需要26.2144ms,這還不包括加法或其它運算，對于大多數實時系統(tǒng)，這個處理時間實在太長。于是尋找一個快速的傅立葉變換算法是人們所期望的。
本來我想把FFT的整個數學推導過程列完出來，但當自己硬著頭皮看完后，發(fā)現對我沒有任何用處，我又不是專門研究數學算法的，哪有那么多時間跟著書本的公式去慢慢推導。我想，這些推導問題還是讓數學家想去吧。我需要的不過是理解它，然后學會應用它就行。有興趣的讀者可以參考相關的資料，這方面的資料實在太多了。
雖然FFT大幅度地降低了常規(guī)傅立葉變換的運算量，但對于一般的單片機而言，處理FFT運算還是力不從心。主要原因是FFT計算過程中的蝶形運算是復數運算，要分開實部和虛部分別計算，想想這是多么繁瑣的事情。可能會有些初學者認為，有這么復雜嗎？我在PC上使用C++一樣可以對復數直接進行加、減、乘、除運算。你說得不錯，可以這么做，但那是C++封裝了對復數處理的類，直接調用就行。在PC上運算這種類型的算法一般不考慮時間和空間，多一兩秒的運行時間不會有什么災難性的結果。
所以我們要衡量一個處理器有沒有足夠的能力來運行FFT算法，根據以上的簡單介紹可以得出以下兩點：
l 處理器要在一個指令周期能完成乘和累加的工作，因為復數運算要多次查表相乘才能實現。其二就是間接尋址，可以實現增/減1個變址量，方便各種查表方法。
l FFT要對原始序列進行反序排列，處理器要有反序間接尋址的能力。
所以，在數字信號的分析處理應用中，DSP比其它的處理器有絕對的優(yōu)勢，因為DSP完全具備以上條件。這就是單片機（51系列，AVR，PIC等等）或ARM處理器很少用來進行數字信號分析的原因。
4． FFT的C實現方法
// 函數名: 快速傅立葉變換（來源《C常用算法集》）
// 本函數測試OK,可以在TC2.0,VC++6.0,Keil C51測試通過。
// 如果你的MCS51系統(tǒng)有足夠的RAM時,可以驗證一下用單片機處理FFT有多么的慢。
//
// 入口參數：
// l: l = 0, 傅立葉變換; l = 1, 逆傅立葉變換
// il: il = 0,不計算傅立葉變換或逆變換模和幅角；il = 1,計算模和幅角
// n: 輸入的點數，為偶數，一般為32，64，128，...,1024等
// k: 滿足n=2^k(k>0),實質上k是n個采樣數據可以分解為偶次冪和奇次冪的次數
// pr[]: l="0時"，存放N點采樣數據的實部
// l="1時", 存放傅立葉變換的N個實部
// pi[]: l="0時"，存放N點采樣數據的虛部
// l="1時", 存放傅立葉變換的N個虛部
//
// 出口參數：
// fr[]: l="0", 返回傅立葉變換的實部
// l="1", 返回逆傅立葉變換的實部
// fi[]: l="0", 返回傅立葉變換的虛部
// l="1", 返回逆傅立葉變換的虛部
// pr[]: il = 1,i = 0 時，返回傅立葉變換的模
// il = 1,i = 1 時，返回逆傅立葉變換的模
// pi[]: il = 1,i = 0 時，返回傅立葉變換的輻角
// il = 1,i = 1 時，返回逆傅立葉變換的輻角
// data: 2005.8.15,Mend Xin Dong
kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0; it<=n-1; it++)
{
m = it;
is = 0;
for(i=0; i<=k-1; i++)
{
j = m/2;
is = 2*is+(m-2*j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
//----------------------------
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306/(1.0*n);
pr[1] = cos(p);
pi[1] = -sin(p);
if (l!=0)
pi[1]=-pi[1];
for (i=2; i<=n-1; i++)
{
p = pr[i-1]*pr[1];
q = pi[i-1]*pi[1];
s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i] = p-q;
pi[i] = s-p-q;
}
for (it=0; it<=n-2; it="it"+2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr+fr[it+1];
fi[it] = vi+fi[it+1];
fr[it+1] = vr-fr[it+1];
fi[it+1] = vi-fi[it+1];
}
m = n/2;
nv = 2;
for (l0=k-2; l0>=0; l0--)
{
m = m/2;
nv = 2*nv;
for(it=0; it<=(m-1)*nv; it="it"+nv)
for (j=0; j<=(nv/2)-1; j++)
{
p = pr[m*j]*fr[it+j+nv/2];
q = pi[m*j]*fi[it+j+nv/2];
s = pr[m*j]+pi[m*j];
s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr = p-q;
poddi = s-p-q;
fr[it+j+nv/2] = fr[it+j]-poddr;
fi[it+j+nv/2] = fi[it+j]-poddi;
fr[it+j] = fr[it+j]+poddr;
fi[it+j] = fi[it+j]+poddi;
}
}
if(l!=0)
{
for(i=0; i<=n-1; i++)
{
fr[i] = fr[i]/(1.0*n);
fi[i] = fi[i]/(1.0*n);
}
}
if(il!=0)
{
for(i=0; i<=n-1; i++)
{
pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if(fabs(fr[i])<0.000001*fabs(fi[i]))
{
if ((fi[i]*fr[i])>0)
pi[i] = 90.0;
else
pi[i] = -90.0;
}
else
pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;
}
}
return;
}
5． 傅立葉變換的幾個重要應用
l 卷積
卷積是濾波網絡對信號響應的術語，即用卷積積分來描述濾波網絡對沖擊函數信號的反應。若x(t)為信號，h(t)為響應，則卷積積分表示如下：
h(t)·x(t) = ∫x(τ)h(t-τ)dτ, 區(qū)間：-∽～+∽
每一個卷積點是信號函數與反轉和平移后的網絡函數的乘積中的區(qū)域。
相關
相關是用于小信號噪聲檢測的一種方法。如果有已知信號與一個噪聲波形相關，用這個方法可以檢測出來，有非零的結果表示發(fā)現了相關性，結果越明顯，相關性越大。
其在形式上與卷積積分相似，如下：
z(t) = ∫x(τ)h(t+τ)dt, 區(qū)間：-∽～+∽
自相關是用來描述一個信號與它自己的相關程度，其值為信號的PSD,即功率譜密度。
1、濾波
這可能是FFT最廣泛的應用了，它使對波形的頻率分量濾波變得十分簡單。比如對采樣信號進行FFT后，干掉不需要的頻率分量，再進行FFT反變換，就得到濾波后的期望信號。
2、信號分析
比如電力監(jiān)控系統(tǒng)的諧波分析，就需要對采樣數據進行FFT運算，然后通過液晶屏或其它人機界面重新繪畫出來，以方便技術人員掌握電力的質量。
小結：
傅立葉變換在目前的相關電子產品中用得非常廣泛，可以說，它是描述函數的另一種語言。掌握傅立葉變換，學會在空域和頻域中同時思考問題，很多時候可以讓我們使用簡單的方法來解決復雜的問題。
