滑動離散傅里葉算法輸出穩定性研究

王宏偉

(西安郵電大學，陜西 西安 710121)

引 言

(1)

式中索引只有頻域索引k。連續計算窗口(窗口長度等于數據幀長度N)滑動后的頻譜，便引入時域索引n，可以進行基于矩形窗的短時傅里葉變換(STFT)的時頻分析，即

n≥N-1

(2)

窗口滑動分為不重疊滑動和重疊滑動。對輸入數據分段(每一段數據稱為一幀數據)的常用方法是不重疊滑動，即從0到N-1，然后從N到2N-1，如此一直下去。另一種方法是對數據幀進行一定程度的重疊處理。例如滑動步長為d時，第一幀取0到N-1點，第二幀取d到d+N-1，第三幀取2d到2d+N-1，依次類推。數據幀最大重疊滑動時為逐點滑動，即第一幀取0到N-1，第二幀取1到N，第三幀取2到N+1，依次類推。

數據逐點滑動時，相鄰窗口內的數據樣本有著很大的相似性，后一個時刻窗口樣本只是將前一個時刻窗口樣本的第一個樣本舍棄，在最后添加一個新的樣本，因此，兩個連續時刻窗口樣本對應的頻譜必然存在聯系，逐點滑動DFT算法[1-4]正是基于這樣的思想進行計算量的優化。

由于數字系統的有限字長效應，存在量化誤差和舍入誤差?；瑒覦FT算法是一種反復迭代的運算，量化誤差和舍入誤差會隨著迭代次數的增加而不衰減的傳播，且隨著時間逐步積累，導致系統失穩。有限字長效應對定點處理數字系統影響的比對浮點處理的數字系統影響更大。

目前解決滑動DFT算法輸出不穩的方法主要有加入衰減因子法[1-2]，自適應最小均方誤差(LMS)算法[3-4]。本文利用改進Goertzel算法的遞歸單元改造滑動DFT算法的頻率單元，使其具有定期自動清零功能，解決了滑動DFT頻率單元輸出不穩定問題。

1 改進Goertzel算法的遞歸單元

圖1(a)為計算DFT譜X(k)的改進Goertzel算法的單極點遞歸單元[5]，其對應的差分方程為

r=1,2,…,N

(3)

式中：r為迭代次數；yk(r)表示第r次迭代結果。令初始狀態yk(0)=0，則

可見，經過N次循環迭代后，得到了我們期望的輸出X(k)=yk(N)。當輸入為復數序列時，每次迭代需要1次復數加法和1次復數乘法，N次迭代共需要N次復數加法和N次復數乘法。

式(3)對應的傳遞函數為

(4)

這是只有一個極點z=ej2πk/N的一階無限脈沖響應(IIR)系統(IIR指具有反饋通道的系統)，但只需進行有限次循環迭代運算。式(4)做如下變換

(ej2πk/N-z-1)

(5)

則Hk(z)=H1(z)H2(z).

H1(z)有兩個極點z=e±j2πk/N，對應的差分方程為

vk(r)= 2cos(2πk/N)·vk(r-1)-

vk(r-2)+x(r-1)

(6)

式中r=1,2,…,N，初始條件為vk(0)=vk(-1)=0.當輸入為實數序列時，式(6)中1次迭代需要1次實數乘法和2次實數加法；N次迭代需要N次實數乘法和2N次實數加法。H2(z)對應差分方程為

yk(r)=vk(r)·ej2πk/N-vk(r-1)

(7)

只需要用式(7)計算當r=N時對應的輸出yk(N)即可

yk(N) =vk(N)·ej2πk/N-v(k)(N-1)

=[vk(N)cos(2πk/N)-vk(N-1)]+

j·vk(N)sin(2πk/N)

=Re[X(k)]+j·Im[X(k)]

=X(k)

(8)

式(8)是在式(6)N次循環迭代完成后再進行實部和虛部的計算，式(8)只需要2次實數乘法和1次實數加法。圖1(b)是實現式(6)和式(8)所示差分方程的雙極點遞歸單元，該遞歸單元計算量為N+2次實數乘法和2N+1次實數加法。由于全部為實數運算，硬件實現結構簡單，計算量少。

圖1所示的兩種改進Goertzel算法的遞歸單元都具有選擇性計算DFT譜X(k)的優點。遞歸單元的輸入序列依次為x(0),x(1),…,x(N-1)，適合于輸入一個數據，便處理一個數據，但要等N次循環迭代后再輸出結果。當輸入為復數數據時，一般用圖1(a)單極點遞歸單元計算DFT譜；而當輸入為實數數據時，選用圖1(b)雙極點遞歸單元計算DFT譜。

(a) 單極點遞歸單元

(b) 雙極點遞歸單元 圖1 改進Goertzel算法的遞歸單元

2 逐點滑動DFT算法介紹

2.1 逐點滑動DFT算法的頻率單元

設在n-1時刻，滑動窗口選中的N個樣本數據為{x(n-N),x(n-N+1),…,x(n-2),x(n-1)}；在n時刻，滑動窗口選中的N個樣本數據為{x(n-N+1),x(n-N+2),…,x(n-1),x(n)}。根據DFT理論，在n-1時刻，第k個頻率單元的頻譜為

X(n-1,k)=x(n-N)+x(n-N+1)e-j2πk/N+

…+x(n-2)e-j2πk·(N-2)/N+

x(n-1)ej2πk·(N-1)/N

(9)

在n時刻，第k個頻率單元的頻譜為

X(n,k)=x(n-N+1)+x(n-N+2)e-j2πk/N+…+

x(n-1)e-j2πk·(N-2)/N+x(n)ej2πk·(N-1)/N

(10)

比較式(9)和(10)可以得到差分方程，即逐點滑動DFT算法的差分方程[1-4]

X(n,k)= [X(n-1,k)-x(n-N)]ej2πk/N+

x(n)e-j2πk(N-1)/N

(11)

式(11)可以化簡為

X(n,k)=ej2πk/N{X(n-1,k)-x(n-N)+x(n)}

(12)

由式(12)可知，要計算n時刻譜值X(n,k)，只要將前一個時刻n-1的譜值X(n-1,k)減去一個舊樣本x(n-N)，加上一個新樣本x(n)，再乘以權系數ej2πk/N即可。圖2(a)是實現式(12)差分方程的單極點頻率單元圖。

式(12)對應的傳遞函數H(z)為

(13)

這是只有一個極點z=ej2πk/N的一階IIR系統(具有反饋通道)，隨著時域索引n→∞，需要進行無限次循環迭代運算。對式(13)再作如下變換：

(14)

vk(n)= 2cos(2πk/N)·vk(n-1)-vk(n-2)+

x(n)-x(n-N)

(15)

H2(z) 對應的差分方程為

X(n,k)=vk(n)·ej2πk/N-vk(n-1)

=[vk(n)cos(2πk/N)-

vK(n-1)]+

j·[vk(n)sin(2πk/N)]

=Re[X(n,k)]+j·Im[X(n,k)]

(16)

圖2(b)是實現式(15)、(16)差分方程的雙極點頻率單元圖。

當輸入序列為復數數據時，一般用圖2(a)單極點頻率單元實現逐點滑動DFT.在一個頻點k處，僅僅需要2次復數加法和1次復數乘法就可以完成頻譜X(n,k)的計算。由于計算機實現時，復數運算需要轉化為實數運算實現，即1次復數與復數的乘法需要4次實數乘法和2次實數加法，1次復數與復數加法需要2次實數加法。因此，圖2(a)單極點頻率單元共需要有4次實數乘法和6次實數加法來完成頻譜X(n,k)的計算。

當輸入為實數數據時，選用圖2(b)雙極點頻率單元實現逐點滑動DFT.在一個頻點k處，僅僅需要3次實數乘法和4次實數加法就可以完成頻譜X(n,k)的計算。由于圖2(b)中全部為實數運算，硬件實現結構簡單，計算量少。

由圖2可知，逐點滑動DFT算法的頻率單元由一個梳狀濾波器(左半部分)和一個復諧振器(右半部分)組合而成[1-2]。如果要同時計算M點頻譜的短時傅里葉變換，只需要M個復諧振器并聯即可。

(a) 單極點頻率單元

(b) 雙極點頻率單元圖2 逐點滑動DFT算法的頻率單元

2.2 逐點滑動DFT算法的優點

1) 逐點滑動DFT算法實現了頻譜輸出的可選擇性。逐點滑動DFT算法采用的是單個頻率單元的獨立譜線計算，根據實際需要，可選擇性計算感興趣的譜線，這極大地提高了譜分析的靈活性。

2) 逐點滑動DFT算法優化了連續對信號進行STFT分析時的計算量。只要知道了前一時刻的頻譜X(n-1,k)，計算下一時刻頻譜X(n,k)時，對于單極點頻率單元，輸入為復數數據時，僅需4次實數乘法、6次實數加法；對于雙極點頻率單元，輸入為實數數據時，僅需3次實數乘法、4次實數加法，并且頻率單元計算量與參數N的大小(即數據幀的長度N)無關。

3) 逐點滑動DFT算法頻率單元的數據輸入率與數據輸出率相同，可以實時處理數據。FFT算法雖然能以較快的速度計算出全部N點離散頻譜，但其N點快速傅里葉變換(FFT)的數據處理時間至少包括N個時域樣本的采集時間和時域數據向頻域數據轉換的計算時間。計算1幀N點FFT時計算量為2Nlog2N次實數乘法和2Nlog2N次實數加法，當數據幀窗口每滑動一個數據后，就需要再獨立的完成1次N點FFT運算，即又需要2Nlog2N次實數乘法和2Nlog2N次實數加法。因此，用FFT算法計算窗口逐點滑動時的STFT運算是不現實的。而對于逐點滑動DFT算法，采集一個時域數據，輸入到頻率單元后，經過幾次算術運算，便可以輸出一個頻域數據。如果在一個采樣間隔時間內，能完成這固定且有限的幾次算術運算，便可以做到邊采樣、邊計算、邊輸出，從而實現了實時處理。

3 遞歸單元和頻率單元的穩定性探討

受數字系統有限字長的影響，表示信號和參加運算的各個參數精度將不再是無限的，存在量化誤差；離散系統中的乘法運算將產生舍入誤差，例如，兩個8 bit的數相乘，其積是16 bit，但最后只能用8 bit來表示，因此，需要舍棄積的后8 bit，這必然產生數字運算的舍入誤差(或截尾誤差)。量化誤差和舍入誤差通過離散系統后，必然產生誤差的積累效應，并在系統中無衰減傳播，結果在系統輸出端將會產生較大的誤差。

(17)

為了導出在遞歸單元輸出端的信噪比，需要對信號作進一步的假設。由N點DFT的定義式(1)，有

(18)

為了防止在運算過程中可能出現的溢出，需要對x(i)的幅度有所限制。不失一般性，假定對所有的i都有|x(i)|≤1，那么|X(k)|≤N；若需要|X(k)|≤1，那么只要將所有的x(i)都除以N即可。為了討論方便，假定x(i)是在(-1/N,1/N)之間服從白噪聲序列，則x(i)的方差為

=1/(3N2)

(19)

于是X(k)的方差為

=1/(3N)

(20)

在遞歸單元輸出端，信號與舍入誤差之間的信噪比為

=6.02b-20lgN

(21)

此處的信噪比(單位dB)指輸出端信號量化表示精度[6]。對于N點DFT的單極點遞歸單元(圖1(a))，若給定DFT的點數N=1 000，并要求SNR=10 dB，由式(21)可以求出所需字長b≈12bit.再如，已知運算單元的字長為20 bit，要求SNR=15 dB，那么由式(21)求出所能允許誤差下的最大迭代次數N≈186 208.

對于N點DFT的雙極點遞歸單元(圖1(b))，求出系數X(k)需要N+2次實數乘法，因此，式(21)需要改寫為

=6.02(b+1)-10lg(N2+2N)

(22)

對于逐點滑動DFT的頻率單元(圖2(a)或圖2(b))，需要進行無限次循環迭代。當n→∞時，相當于式(21)或(22)中的N趨于無窮大，在有限字長運算單元組成的數字系統中，隨著循環迭代次數的無限增加，信噪比越來越小，意味著輸出端信號值的精度越來越差，導致輸出結果不可靠。由于量化誤差和舍入誤差的積累與傳播，逐點滑動DFT頻率單元輸出不穩定現象遲早會出現。

4 改進滑動DFT算法的頻率單元

由上節討論可知，改進Goertzel算法的遞歸單元和逐點滑動DFT算法的頻率單元雖然都采用了具有反饋通路的IIR濾波器結構，但遞歸單元進行有限次循環迭代，在大字長運算單元組成的遞歸單元中，確實可以保證輸出穩定；而頻率單元則需要進行無限次循環迭代運算，因長時間的誤差積累會出現輸出不穩的致命缺點。

利用改進Goertzel算法的單極點遞歸單元(圖1(a))來改造逐點滑動DFT算法的單極點頻率單元(圖2(a))，改進滑動DFT算法的單極點頻率單元，如圖3(a)所示。圖3(a)上半部分為改進Goertzel算法的單極點遞歸單元，每輸入N個點才會有一次輸出Xr(lN-1,k)；下半部分為逐點滑動DFT算法的單極點頻率單元，每輸入一個點便有一次輸出Xs(n,k)。圖3(a)中選擇器模塊的輸出滿足

X(n,k)=

(23)

利用改進Goertzel算法的雙極點遞歸單元(圖1(b))來改造逐點滑動DFT算法的雙極點頻率單元(圖2(b))，改進滑動DFT算法的雙極點頻率單元，如圖3(b)所示。圖3(b)上半部分為改進Goertzel算法的雙極點遞歸單元的反饋部分，每輸入N個點才會有輸出vr(lN-1,k)，vr(lN-2,k) ；下半部分為逐點滑動DFT算法的雙極點頻率單元，每輸入一個點便有輸出vs(n,k)，vs(n-1,k)。圖3(b)中選擇器1模塊的輸出滿足

(a) 改進滑動DFT的單極點頻率單元

(b) 改進滑動DFT雙極點頻率單元 圖3 改進滑動DFT頻率單元

vk(n)=

(24)

選擇器2模塊的輸出滿足

(25)

由于單極點遞歸單元(圖1(a))在N次循環迭代運算前需要令yk(0)=0；雙極點遞歸單元(圖1(b))在循環迭代運算前也需要令vk(0)=vk(-1)=0.計算機具體實現yk(0)=0或vk(0)=v(-1)=0時，需要對乘法器、加法器和鎖存器等運算單元執行清零操作。改進Goertzel遞歸單元的自動清零操作防止了量化誤差和舍入誤差的進一步積累和傳播。

由于遞歸單元有定期自動清零功能，則滑動DFT的頻率單元的輸出每隔N個值，便被遞歸單元的輸出值替代一次，相當于定期給滑動DFT頻率單元清零，并能提供準確的新譜值，所以改進滑動DFT的頻率單元(圖3)防止了誤差的積累和傳播。只要保證N次循環迭代運算中，遞歸單元的輸出穩定(采用大字長的定點或浮點運算單元，擴大數值表示范圍，可以保證輸出穩定)，則可保證改進滑動DFT頻率單元的輸出穩定。

5 計算機仿真

為了方便對比分析，以文獻[4]列舉的信號為例研究。信號x(n)=sin(2πF1·n/fs)+sin(2πF2·n/fs)，其中F1=0.03 Hz；F2=1.1 Hz；采樣頻率fs=32 Hz.

5.1 首幀數據對應的頻譜圖計算

比較FFT算法和改進Goertzel算法得到的首幀數據對應的頻譜圖(圖4)，無論是幅頻圖，還是相頻圖都完全相同，從而驗證了這兩種算法都可以被用來計算信號的DFT譜。但FFT算法不管需要與否，會不加選擇的計算出全部離散頻譜，適合全景頻譜分析；而改進Goertzel算法采取單根譜線的獨立計算，具有選擇性計算頻譜的優點，適合感興趣頻點的頻譜計算。

不論輸入數據是實信號，還是復信號，N點FFT算法的計算量均指復數乘法和復數加法，那么32點FFT需要(N/2)log2N=16log232=80次復數乘法和Nlog2N=32log232=160次復數加法。

圖4 兩種算法得到的頻譜圖比較

輸入實信號時，如果以1路改進Goertzel算法的雙極點遞歸單元(圖1(b))為運算模塊，依次計算全部N=32點頻譜時，計算量為N(N+2)=1 026次實數乘法和N(2N+1)=2 080次實數加法。當然，若考慮了某些頻點的特殊性和DFT譜的對稱性還可以進一步減少計算量。例如根據式(1)，有

X(0) =X(k)∣k=0

(26)

(27)

=x(0)-x(1)+x(2)-x(3)+…+

(-1)N-1x(N-1)

本文將改進Goertzel算法的遞歸單元與滑動DFT算法的頻率單元集成在一起構成改進滑動DFT算法的頻率單元，只負責一個指定頻點的DFT譜計算；如果需要計算多個指定頻點或局部頻段的多個頻點DFT譜，一般需要多路頻率單元構成并行運算模塊。此時模塊的計算量(或計算速度)只需考慮單路遞歸單元的計算量(或計算速度)即可。對于本例，單路遞歸單元計算量為N+2=34次實數乘法和2N+1=65次實數加法。與FFT算法不同，改進Goertzel算法不必等參與運算的數據都準備好，可以邊采集數據、邊計算頻譜，具有很好的實時性。

5.2 有限字長效應對滑動DFT輸出影響

作N=32點的逐點滑動DFT。將每個時刻m對應DFT譜各個分量的模平方(功率)累加起來，除以N作為輸出，即

(28)

Matlab軟件是一種64位浮點計算引擎，并給出十進制數的結果[9]，一般短期內不會因為有限字長效應出現輸出不穩定現象。而現有數字信號處理(DSP)器件或現場可編程門陣列(FPGA)器件上采用有限字長(字長為8 bit,12 bit等)的定點或浮點二進制運算則可能出現輸出不穩定現象，因此，需要用Matlab仿真有限字長的定點或浮點二進制運算。仿真時需要對輸入數據、各種運算參數及中間結果均做量化處理，并對乘法運算結果做量化舍入處理后，研究有限字長對滑動DFT算法輸出穩定性的影響。本文仿真時取字長b=8 bit的定點運算單元。

仿真時，將不作量化和舍入處理的滑動DFT算法的輸出結果作為理想輸出，與經過量化和舍入處理的滑動DFT算法輸出結果，以及與經過量化和舍入處理并經過LMS算法(文獻[4]中介紹的算法)改造后的滑動DFT算法輸出結果，還有與經過量化和舍入處理并經過改進Goertzel算法改造后的滑動DFT算法(圖3(b)頻率單元)輸出結果，進行比較，如圖5所示，圖中橫坐標代表時間序列m.

圖5 經過不同處理后的滑動DFT輸出結果比較

在圖5中，經過量化和舍入處理的滑動DFT算法輸出結果隨著時間的推移，出現了"漂移現象"(誤差積累導致失穩)。而經過量化和舍入處理并經過LMS算法改造后的滑動DFT算法輸出結果，與經過量化和舍入處理并經過改進Goertzel算法改造后的滑動DFT算法輸出結果，都不會出現"漂移現象"，可以長時間連續分析信號時頻譜。但經過改進Goertzel改造后的滑動DFT算法，與經過LMS算法改造后的滑動DFT算法比較，計算量更少，控制邏輯更簡單。

6 結 論

由改進Goertzel算法的遞歸單元對逐點滑動DFT算法的頻率單元改造后，改進滑動DFT算法的頻率單元不會出現輸出不穩定現象，這種方法在連續實時頻譜分析中具有重要的意義。

由于時域加窗(指非矩形窗)會破壞數據幀滑動的連續性，致使滑動DFT算法的時域差分方程不再成立，需要在頻域加窗，以減少頻譜泄露。關于頻域加窗代替時域加窗，以減少頻譜泄露方面的內容，可以參考文獻[1-4]，本文不再贅述。

[1] JACOBSEN E, LYONS R. The sliding DFT[J] . IEEE Signal Processing Magazine, 2003, 20(3): 74-80.

[2] JACOBSEN E, LYONS R. An update to the sliding DFT[J]. IEEE Signal Processing Magazine, 2004, 21(1): 110-111.

[3] FARHANG-BOROUJENY B, GAZOR S. Generalized sliding FFT and its application to implementation of block LMS adaptive filters[J]. IEEE Transactions on Signal Processing, 1994, 42(3): 532-538.

[4] BEAUFAYS F, WIDROW B. On the advantages of the LMS spectrum analyzer over the nonadaptive implementations of the sliding-DFT[J]. IEEE Transaction on Circuits and Systems I: Fundamental Theory and Applications, 1995, 42(4): 218-220.

[5] 王宏偉, 趙國慶. 遞歸算法的參數設置[J].電波科學學報, 2010, 25(6): 1187-1190.

WANG Hongwei, ZHAO Guoqing. Parameter setting of recursive algorithm[J]. Chinese Journal of Radio Science, 2010, 25(6): 1187-1190. (in Chinese)

[6] 胡廣書. 數字信號處理——理論、算法與實現[M]. 北京: 清華大學出版社, 2003.

[7] GERALD G.An algorithm for the evaluation of finite trigonometric series[J]. The American Mathematical Monthly, 1958, 65(1): 34-35.

[8] BANKS K. The Goertzel algorithm[J]. Embedded System Programming, 2002, 15(9): 32-42.

[9] 維納 K 恩格爾, 約翰G普羅克斯. 數字信號處理(使用MATLAB)[M]. 劉樹棠, 譯. 西安: 西安交通大學出版社, 2002.