郎銳 頻譜分析是電子工程上一個非常重要的分析手段，許多計算機(jī)輔助電路分析（CAA）類軟件都具備這種分析能力，以便電子工程師能清楚地看到某波形的頻譜分布情況，
用VC＋＋實現(xiàn)對波形數(shù)據(jù)的頻譜分析.net
。要對一個輸入信號源作頻譜分析，將其由時域信號轉(zhuǎn)變?yōu)轭l域信號，就必然要用到傅立葉變換。 這
郎銳
頻譜分析是電子工程上一個非常重要的分析手段，許多計算機(jī)輔助電路分析（CAA）類軟件都具備這種分析能力，以便電子工程師能清楚地看到某波形的頻譜分布情況。要對一個輸入信號源作頻譜分析，將其由時域信號轉(zhuǎn)變?yōu)轭l域信號，就必然要用到傅立葉變換。這樣，無論是在時域還是在頻域，都要對連續(xù)函數(shù)進(jìn)行積分運算。很顯然，要通過計算機(jī)實現(xiàn)這種變換就需要預(yù)先通過抽樣將原始的連續(xù)數(shù)據(jù)轉(zhuǎn)變?yōu)殡x散數(shù)據(jù)，并將計算范圍收縮到一個有限區(qū)間。因此，在允許一定程度近似的條件下，可以使用“離散傅立葉變換（DFT）”對波形數(shù)據(jù)進(jìn)行頻譜分析。
算法構(gòu)成原理
要計算一個N點的離散傅立葉變換需要同一個N×N點的W矩陣（關(guān)于W矩陣請參閱信號與系統(tǒng)方面或數(shù)學(xué)方面的書籍）相運算，隨著N值的增大，運算次數(shù)顯著上升，當(dāng)點數(shù)達(dá)到1024時，需要進(jìn)行復(fù)數(shù)乘法運算1048576次。顯然這種算法在實際運用中無法保證當(dāng)點數(shù)較大時的運算速度，無法滿足對信號的實時處理要求。
根據(jù)W矩陣中W元素的周期性和對稱性我們可以將一個N點的DFT運算分解為兩組N/2點的DFT運算，然后取和即可。為進(jìn)一步提高效率，將上述兩個矩陣按奇偶順序逐級分解下去。當(dāng)采樣點數(shù)為2的指數(shù)次方M時，可分解為M級子矩陣運算，全部工作量如下：
復(fù)數(shù)乘法：M×N/2次
復(fù)數(shù)加法：N×M次
直接采用DFT算法需要的運算量為：
復(fù)數(shù)乘法：N×N次
復(fù)數(shù)加法：N×(N－1)次
當(dāng)點數(shù)N為幾十個點時快速傅立葉交換(FFT)的優(yōu)勢還不明顯，而一旦N達(dá)到幾千時優(yōu)勢是十分明顯的：
N=1024時：DFT需1048576次運算，F(xiàn)FT僅需5120次運算，改善比為204.8。
N=2048時：DFT需4194304次運算，F(xiàn)FT僅需11264次運算，改善比達(dá)到372.4。
當(dāng)采樣點數(shù)較多時，如變換前和變換后的序列都按自然順序排列，則中間運算過程會占用大量的中間存儲單元，造成效率的低下和存儲單元的浪費。根據(jù)FFT的實現(xiàn)原理我們可以對采樣序列進(jìn)行逐次奇偶抽選，打亂以前的次序重新排序，然后按此順序參加運算，以“即位運算”提高存儲單元的利用率。
復(fù)數(shù)的描述方法
進(jìn)行傅立葉變換時不可避免地要用到復(fù)數(shù)，而在VC中并沒有現(xiàn)成的可用于表示復(fù)數(shù)的數(shù)據(jù)類型，因此需要自己定義一個含有兩個成員變量的數(shù)據(jù)結(jié)構(gòu)來表示復(fù)數(shù)，這兩個成員變量可分別用于表示復(fù)數(shù)的實部與虛部：
typedef struct tagComplex{
//復(fù)數(shù)的實部
float Re;
//復(fù)數(shù)的虛部
float Im;
}Complex;
倒序的實現(xiàn)
在進(jìn)行快速傅立葉變換時，可以將輸入的時域序列和輸出的頻域序列都按照自然順序排列；也可以按照“蝴蝶圖”所描述的計算方法對輸入的時域序列按奇偶分解后的序列排序，而輸出的頻域序列仍是按自然順序排列；還有一種方式是輸入的時域序列是自然序列，而輸出的頻域序列則是按奇偶分解后的順序排列。這三種方式各有優(yōu)缺點：第一種對輸入、輸出不需要進(jìn)一步排序，但由于自然排序不符合“蝴蝶圖”運算規(guī)律，會占用大量中間存儲單元；而后兩種則無需中間存儲單元，但需要倒序。權(quán)衡利弊，當(dāng)采樣點較多時還是采用后兩種方式好，多一次倒序運算對現(xiàn)在的高性能計算機(jī)而言并不是什么負(fù)擔(dān)。下面代碼用于對原始采樣序列的時間抽選奇偶分解工作，其中A、N分別表示指向采樣序列復(fù)數(shù)數(shù)組的指針和序列的長度。
int NV2=N/2;
int NM1=N－1;
int I,J,K=0;
//用于中介的復(fù)數(shù)變量T
Complex T;
I=J=1;
while(I<=NM1)
{
if(I< J)
{
//借助于中間變量T，將A[J－1]的內(nèi)容和A[I－1]的內(nèi)容互換
T=A[J－1];
A[J－1]=A[I－1];
A[I－1]=T;
}
K=NV2;
while(K< J)
{
J－=K;
K/=2;
}
J＋=K;
I＋＋;
}
時域信號的頻譜分析
首先要將從外設(shè)輸入或采集的時域波形數(shù)據(jù)經(jīng)抽樣量化后，通過CFile類的Open（……）、Read(……)等成員函數(shù)將其讀取到緩存中，并將其轉(zhuǎn)化為復(fù)變量存放于復(fù)變量數(shù)組A中。同時需要驗證數(shù)據(jù)量的長度是否為2的整數(shù)次冪，如不是則用0來補齊，否則無法用“蝴蝶圖”進(jìn)行分解運算。下面代碼用于完成對原始采樣時域序列的快速傅立葉變換，A、M分別表示指向原始采樣數(shù)據(jù)數(shù)組的指針和序列長度的2的整數(shù)次冪：
……
Complex U,W,T;
int LE,LE1,I,J,IP;
int N=(int)pow(2,M);
//由于采用時間抽選奇偶分解方式，所以在參加運算前首先要對時間序列進(jìn)行倒序
ReverseOrder(A,N);
int L=1;
while(L<=M)
{
LE=(int)pow(2,L);
LE1=LE/2;
U.Re=1.0f;
U.Im=0.0f;
//計算W算子的值
W.Re=(float)cos(PI/(1.0＊LE1));
W.Im=(float)－1.0＊sin(PI/(1.0＊LE1));
if(abs(W.Re)<1.0e－12)
W.Re=0.0f;
if(abs(W.Im)< 1.0e－12)
W.Im=0.0f;
J=1;
while(J<=LE1)
{
I=J;
while(I<=N)
{
IP=I＋LE1;
//復(fù)數(shù)運算A×U
T.Re=(float)A[IP－1].Re＊U.Re－A
[IP－1].Im＊U.Im;
T.Im=(float)A[IP－1].Re＊U.Im＋A
[IP－1].Im＊U.Re;
//復(fù)數(shù)運算A－T
A[IP－1].Re=(float)A[I－1].Re－
T.Re;
A[IP－1].Im=(float)A[I－1].Im－
T.Im;
//復(fù)數(shù)運算A＋T
A[I－1].Re＋=T.Re;
A[I－1].Im＋=T.Im;
I＋=LE;
}
float temp=U.Re;
//復(fù)數(shù)運算U×W
U.Re=(float)U.Re＊W.Re－U.Im＊
W.Im;
U.Im=(float)temp＊W.Im＋U.Im＊
W.Re;
J＋＋;
}
L＋＋;
}
……
上述代碼執(zhí)行完畢時，原先存放著時域數(shù)值的復(fù)變量數(shù)組內(nèi)存放的就是經(jīng)過分析后的頻域值，利用此數(shù)據(jù)可以通過繪圖將頻域波形直觀地顯示出來，也可以將其存成數(shù)據(jù)文件，以備進(jìn)一步使用，
電腦資料
《用VC＋＋實現(xiàn)對波形數(shù)據(jù)的頻譜分析.net》(http://www.lotusphilosophies.com)。測試及運算結(jié)果分析
編譯運行程序，分析一個三角脈沖的數(shù)據(jù)文件，并保存分析結(jié)果。該三角脈沖幅度為1，持續(xù)時間2毫秒，抽樣時間間隔是20微秒，延拓周期（數(shù)據(jù)記錄長度）為10毫秒，采樣點數(shù)為500，取2的整數(shù)次冪512個采樣點。下附該三角脈沖頻譜的計算結(jié)果及誤差分析：
頻率（Hz） FFT求得 X(f) 誤差
0.00 1.00006E－03 1.00000E－03 6.10352E－08
100.00 9.67593E－04 9.67531E－04 6.14332E－08
200.00 8.75203E－04 8.75150E－04 6.25092E－08
300.00 7.36904E－04 7.36849E－04 6.39413E－08
400.00 5.72852E－04 5.72787E－04 6.52926E－08
500.00 4.05351E－04 4.05285E－04 6.61362E－08
600.00 2.54638E－04 2.54572E－04 6.61847E－08
……
2700.00 9.16539E－06 9.09679E－06 6.86075E－08
2800.00 4.53216E－06 4.46500E－06 6.71550E－08
2900.00 1.21487E－06 1.15945E－06 6.44190E－08
注：此處FFT運算結(jié)果都乘以了系數(shù)10毫秒（0.01秒）。
從上述數(shù)據(jù)中可以看出，在分析結(jié)果中產(chǎn)生了誤差。這是由于待分析的連續(xù)時間信號不具備離散性或周期性，也可能有無限長度。為了適應(yīng)FFT方法的需要，先對波形進(jìn)行了抽樣和截斷，這樣再用程序分析采樣數(shù)據(jù)必然會引起誤差。從分析結(jié)果還可以看出，頻率越高，誤差波動也越大，此分析結(jié)果產(chǎn)生的誤差在允許范圍之內(nèi)，是一個可以允許的近似。
本程序在Windows98、Microsoft Visual C＋＋ 6.0下編譯通過。
原文轉(zhuǎn)自：http://www.ltesting.net