用VC＋＋實現對波形數據的頻譜分析

郎銳

頻譜分析是電子工程上一個非常重要的分析手段，許多計算機輔助電路分析（CAA）類軟件都具備這種分析能力，以便電子工程師能清楚地看到某波形的頻譜分布情況。要對一個輸入信號源作頻譜分析，將其由時域信號轉變為頻域信號，就必然要用到傅立葉變換。這樣，無論是在時域還是在頻域，都要對連續函數進行積分運算。很顯然，要通過計算機實現這種變換就需要預先通過抽樣將原始的連續數據轉變為離散數據，并將計算范圍收縮到一個有限區間。因此，在允許一定程度近似的條件下，可以使用“離散傅立葉變換（DFT）”對波形數據進行頻譜分析。  
算法構成原理
要計算一個N點的離散傅立葉變換需要同一個N×N點的W矩陣（關于W矩陣請參閱信號與系統方面或數學方面的書籍）相運算，隨著N值的增大，運算次數顯著上升，當點數達到1024時，需要進行復數乘法運算1048576次。顯然這種算法在實際運用中無法保證當點數較大時的運算速度，無法滿足對信號的實時處理要求。
根據W矩陣中W元素的周期性和對稱性我們可以將一個N點的DFT運算分解為兩組N/2點的DFT運算，然后取和即可。為進一步提高效率，將上述兩個矩陣按奇偶順序逐級分解下去。當采樣點數為2的指數次方M時，可分解為M級子矩陣運算，全部工作量如下：
復數乘法：M×N/2次
復數加法：N×M次
直接采用DFT算法需要的運算量為：
復數乘法：N×N次
復數加法：N×(N－1)次
當點數N為幾十個點時快速傅立葉交換(FFT)的優勢還不明顯，而一旦N達到幾千時優勢是十分明顯的：
N=1024時：DFT需1048576次運算，FFT僅需5120次運算，改善比為204.8。
N=2048時：DFT需4194304次運算，FFT僅需11264次運算，改善比達到372.4。
當采樣點數較多時，如變換前和變換后的序列都按自然順序排列，則中間運算過程會占用大量的中間存儲單元，造成效率的低下和存儲單元的浪費。根據FFT的實現原理我們可以對采樣序列進行逐次奇偶抽選，打亂以前的次序重新排序，然后按此順序參加運算，以“即位運算”提高存儲單元的利用率。
復數的描述方法
進行傅立葉變換時不可避免地要用到復數，而在VC中并沒有現成的可用于表示復數的數據類型，因此需要自己定義一個含有兩個成員變量的數據結構來表示復數，這兩個成員變量可分別用于表示復數的實部與虛部：
typedef struct tagComplex{
//復數的實部
float Re;  
//復數的虛部
float Im;  
}Complex;
倒序的實現
在進行快速傅立葉變換時，可以將輸入的時域序列和輸出的頻域序列都按照自然順序排列；也可以按照“蝴蝶圖”所描述的計算方法對輸入的時域序列按奇偶分解后的序列排序，而輸出的頻域序列仍是按自然順序排列；還有一種方式是輸入的時域序列是自然序列，而輸出的頻域序列則是按奇偶分解后的順序排列。這三種方式各有優缺點：第一種對輸入、輸出不需要進一步排序，但由于自然排序不符合“蝴蝶圖”運算規律，會占用大量中間存儲單元；而后兩種則無需中間存儲單元，但需要倒序。權衡利弊，當采樣點較多時還是采用后兩種方式好，多一次倒序運算對現在的高性能計算機而言并不是什么負擔。下面代碼用于對原始采樣序列的時間抽選奇偶分解工作，其中A、N分別表示指向采樣序列復數數組的指針和序列的長度。
int NV2=N/2;
int NM1=N－1;
int I,J,K=0;
//用于中介的復數變量T
Complex T;
I=J=1;
while(I<=NM1)
{
if(I< J)
{
//借助于中間變量T，將A[J－1]的內容和A[I－1]的內容互換
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＋＋;
}
時域信號的頻譜分析
首先要將從外設輸入或采集的時域波形數據經抽樣量化后，通過CFile類的Open（……）、Read(……)等成員函數將其讀取到緩存中，并將其轉化為復變量存放于復變量數組A中。同時需要驗證數據量的長度是否為2的整數次冪，如不是則用0來補齊，否則無法用“蝴蝶圖”進行分解運算。下面代碼用于完成對原始采樣時域序列的快速傅立葉變換，A、M分別表示指向原始采樣數據數組的指針和序列長度的2的整數次冪：
……
Complex U,W,T;
int LE,LE1,I,J,IP;
int N=(int)pow(2,M);
//由于采用時間抽選奇偶分解方式，所以在參加運算前首先要對時間序列進行倒序
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;
//復數運算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;
//復數運算A－T
A[IP－1].Re=(float)A[I－1].Re－
T.Re;
A[IP－1].Im=(float)A[I－1].Im－
T.Im;
//復數運算A＋T
A[I－1].Re＋=T.Re;
A[I－1].Im＋=T.Im;
I＋=LE;
}
float temp=U.Re;
//復數運算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＋＋;
}
……
上述代碼執行完畢時，原先存放著時域數值的復變量數組內存放的就是經過分析后的頻域值，利用此數據可以通過繪圖將頻域波形直觀地顯示出來，也可以將其存成數據文件，以備進一步使用。
測試及運算結果分析
編譯運行程序，分析一個三角脈沖的數據文件，并保存分析結果。該三角脈沖幅度為1，持續時間2毫秒，抽樣時間間隔是20微秒，延拓周期（數據記錄長度）為10毫秒，采樣點數為500，取2的整數次冪512個采樣點。下附該三角脈沖頻譜的計算結果及誤差分析：
頻率（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運算結果都乘以了系數10毫秒（0.01秒）。
從上述數據中可以看出，在分析結果中產生了誤差。這是由于待分析的連續時間信號不具備離散性或周期性，也可能有無限長度。為了適應FFT方法的需要，先對波形進行了抽樣和截斷，這樣再用程序分析采樣數據必然會引起誤差。從分析結果還可以看出，頻率越高，誤差波動也越大，此分析結果產生的誤差在允許范圍之內，是一個可以允許的近似。
本程序在Windows 98、Microsoft Visual C＋＋ 6.0下編譯通過。 

原文轉自：http://www.anti-gravitydesign.com