基于嵌入式系統(tǒng)的頻率計權實現(xiàn)方法分析與對比

羅自榮，陳章位，蔡德威，樊亞容

(浙江大學 流體傳動與機電系統(tǒng)國家重點實驗室,杭州 310013)

近年來人們對噪聲污染的關注不斷加大,聲學分析技術在噪聲評估,故障診斷,污染防治,消噪減噪,疾病治療等方面的應用也不斷加深[1-2]。人耳對聲音的感受不僅跟聲壓大小有關,而且跟聲壓的頻率有關。在聲學分析中使用頻率計權網(wǎng)絡,以達到模擬人耳聽覺的效果。

模擬電路時代,頻率計權的實現(xiàn)主要是通過利用電容,電阻,運放電路等組成計權網(wǎng)絡實現(xiàn)頻率計權,但是這種方法由于各種電子器件的精度較難保證,而且設計比較復雜,不易變更和調試,同時可靠性和穩(wěn)定性都比較低。近年來全數(shù)字式的頻率計權方法得到廣泛的應用[3],通過設計一個頻率計權濾波器,然后把采樣后的聲壓信號與濾波器做卷積運算,得到計權結果,但隨著濾波精度的提高,濾波器系數(shù)的增長,計權運算需要的運算量也急劇增大。本文提出了基于重疊相加DFT的優(yōu)化方案,并對實現(xiàn)過程進行理論分析。對基于卷積和基于重疊相加DFT的方案的運算量級進行理論推導,計算和討論。理論推導和實驗結果表明：兩種方案的結果是一致的,隨著濾波器系數(shù)增大后者的性能快速提升,并得出在濾波器系數(shù)已知,和濾波器系數(shù),一幀大小都已知的情況下的方案和參數(shù)選擇準則。

1 基于卷積的全數(shù)字式頻率計權實現(xiàn)

1.1 頻率計權聲壓

在聲學測量儀器中,依據(jù)不同需求使用不同的計權網(wǎng)絡,現(xiàn)在常用的頻率計權網(wǎng)絡有A,B,C 計權。

圖1 A,B,C計權網(wǎng)絡衰減曲線

圖1所示的為A,B和C計權網(wǎng)絡衰減曲線圖?？梢钥吹饺硕鷮Ω哳l段和低頻段的聲壓的反應比較遲鈍,對1 000 Hz左右的聲壓反應比較靈敏。式(1)和(2)分別是A,C計權頻率與對應的衰減量的函數(shù)關系式。其中：

f1=20.6 Hz,f2=107.7 Hz,

f3=737.9 Hz,f4=12 194 Hz,

A1 000=2.000 dB,C1 000=-0.063 dB

A(f)=

1.2 基于卷積的頻率計權實現(xiàn)

聲壓頻率計權可通過讓聲壓信號經(jīng)過對應頻率計權濾波器實現(xiàn)。使用雙線性變換法,結合(1)式得A,C計權所對應的頻率計權濾波器的傳遞函數(shù)[4]

(2)

對頻率計權濾波器傳遞函數(shù)從S平面映射到Z平面得[5]：

(3)

其中:a0=1,a1=-0.984 6,a2=-0.012 65,b0=0.505 7,b1=-0.505 7,a10=1,a11=-1.894 4,a12=0.895 7,b10=0.947 5,b11=-1.895,b12=0.974 5。

得到濾波器傳遞函數(shù)Z變換,容易得濾波器的系數(shù),假設取濾波器的系數(shù)長度為M,調用卷積計算式(3)即可以求得對應頻率計權的結果[6]。

y(n)=x(n)*h(n)

(4)

如此一幀一幀處理,便可得頻率計權后的信號。但是隨著分析精度要求的提高,逐漸要求濾波器系數(shù)越來越長,而計權濾波的運算量跟濾波器長度有關。容易得到對于一列長為L的數(shù)據(jù),通過一個系數(shù)長度為M的濾波器總共需要L×M次乘除法和L×(M-1)次加減法運算。所以當M增大時,處理需要的運算量也急劇的增大。

2 基于重疊相加DFT的頻率計權實現(xiàn)

2.1 基于重疊相加DFT的卷積實現(xiàn)

隨著計權濾波器系數(shù)長度的增大,運算量急劇上升,考慮用循環(huán)卷積計算線性卷積,而依據(jù)卷積定理可通過DFT/IDFT計算循環(huán)卷積[7-9]。設濾波器系數(shù)h(n)長度為M,一幀長度為N的數(shù)據(jù)x(n)(x(n)是N個時域聲壓信號以N為周期的周期延拓,是確定信號,可以使用DFT進行分析),yc(n)和y(n)分別表示h(n)與x(n)的L點循環(huán)卷積和線性卷積,且滿足L>max[N,M],則有:

yc(n)=h(n)?x(n)=

(5)

x((n-m))L為x(n)的以L為周期的周期延拓,RL(n)為單位矩形序列,當0≤n≤L-1時RL(n)=1,其他情況RL(n)=0所以有：

(6)

又依卷積的定義有：

(7)

聯(lián)立式(5),式(6)和式(7)得：

(8)

所以h(n)和x(n)的L點循環(huán)卷積為線性卷積y(n)以L為周期的周期延拓序列的主值序列。由定義可知y(n)序列的長度為N+M-1。所以y(n)以L為周期的周期延拓無混疊的條件是:

L≥N+M-1

(9)

即當循環(huán)卷積長度L,一幀長度N和濾波器系數(shù)長度M滿足式(5)時有：yc(n)=y(n)。此時依卷積定理通過DFT/IDFT計算得h(n)和x(n)的L點循環(huán)卷積,進而得到線性卷積值。具體計算過程如圖2所示,先對h(n)和x(n)分別補L-N和L-M個零,得L點序列。然后分別進行L點DFT得H(k)和X(k)。H(k)乘X(k)得Yc(k)。對Yc(k)進行L點IFFT得yc(n)。yc(n)便是h(n)和x(n)的線性卷積值y(n)。

圖2 卷積計算流程圖

這里采用重疊相加DFT計算得h(n)和x(n)的卷積,避免卷積計算過程中繁瑣重復的相乘相加運算,而對DFT,IDFT可以采用快速算法FFT和IFFT計算求得。

對于無限長序列x(n),與先前類似,分成一幀一幀,每幀長度N,設xi(n)=x(n+iN)RN(n),則：

(10)

設yi(n-iN)=h(n)*xi(n-iN),有：

y(n)=h(n)*x(n)=

(11)

從式(11)可知對于無限長序列x(n)經(jīng)過頻率計權后的結果y(n)為各幀卷積重疊相加。當L=N+M-1時,即為第(i-1)幀的末尾M-1個值與第i幀的前面(M-1)個值相加。具體過程如圖3所示。

圖3 重疊相加法計算卷積步驟

2.2 頻率計權中DFT的實現(xiàn)

采用FFT計算DFT很大得減少了DFT的運算量,但是FFT運算是針對復數(shù)的,而要處理的聲壓信號都是實數(shù),實數(shù)可以看為虛部為0的復數(shù),但很明顯這有很多的冗余運算,所以考慮把兩個相鄰的實數(shù)看為一個復數(shù)[10],對于L點序列x(n),設x1(n)和x2(n)分別為序列的奇序列和偶序列即g(n)=x1(2n)+jx2(2n+1)。所以G(k)=X1(k)+jX2(k),則有:

(12)

又結合DFT的定義和性質有：

(13)

式(12)與式(13)聯(lián)立有：

(14)

所以要計算x(n)的L點DFT,只需計算g(n)的L/2點DFT然后利用式(14)即可計算得X(k)的值。

2.3 頻率計權中IDFT的實現(xiàn)

按圖2計算得Yc(k),需要進行L點IDFT才能求得最終的頻率計權結果。類似DFT的運算,是否有一種簡便的方法減少IDFT需要計算的點數(shù)?？紤]頻率計權的聲壓信號是實數(shù),而設計的濾波器的系數(shù)也都是實數(shù),實數(shù)與實數(shù)的卷積的結果肯定也是實數(shù)。也就是說在對Yc(k)進行L點IDFT后的結果是為L個實數(shù)。即Yc(k)為yc(n)的L點實序列DFT,設Y1(k)和Y2(k)分別為yc(n)的奇序列和偶序列DFT。則有[11]：

(15)

(16)

另一方面，yc(n)的派生復數(shù)序列,即把相鄰的兩個實數(shù),一個作為復數(shù)的實部,一個作為復數(shù)的虛部，g(n)=y1(n)+jy2(n),即G(k)=Y1(k)+jY2(k),則聯(lián)立式(15)和式(16)有：

G(k)=Y1(k)+jY2(k)=

(17)

所以在頻率計權時,對于一L點DFTYc(k)要求它的IDFT,只需要按照式(17)求得G(k),然后對G(k)做L/2點IDFT得g(n),g(n)的實部就是yc(n)的奇序列,虛部就是yc(n)的偶序列。

3 算法復雜度分析與對比

設嵌入式系統(tǒng)做一個乘除法運算的時間為Tm,一個加減法運算需要的時間為Ta。濾波器的系數(shù)的長度為M,考慮嵌入式系統(tǒng)做卷積實現(xiàn)的需要和避免頻繁的數(shù)據(jù)移動,一般一幀的長度N-1應大于或等于M,則對于直接采用卷積定義的頻率計權實現(xiàn)方法計算一列長為Le(Le→∞)的序列x(n)需要的運算量為

(18)

對于基于重疊相加DTF的優(yōu)化方案,設一幀的長度為N,Le=nN,取循環(huán)卷積長度L=N+M-1,則需要的運算量為

(19)

在DSP中采用了硬件乘法器的結構,乘法運算的運算量大大減小,跟加法指令的運算周期接近,這里取Ta=Tm=T,所以：

(20)

對式(20)分子分母同時除以M得：

對一幀大小已知的情況,即N為已知量,m為關于M的函數(shù),顯然式(21)的分子為關于變量M的單調增函數(shù)。對于式(21)的分母,取k=N-1,則分母可表示為以M為自變量的函數(shù)f(x)：

對f(x)求導得：

(22)

考慮前面分析N-1≥M,即x≤k,且濾波器系數(shù)M大于1,即：x>1,所以對于式(22)有：

所以f(x)為一個單調減函數(shù),對于式(21),分子是M的單調增函數(shù),分母是M的單調減函數(shù),所以M越大,基于重疊相加DFT優(yōu)化算法的優(yōu)越性越好。

對于式(20),分子分母同時除以N得:

(23)

對于濾波器設計已經(jīng)完成,即M已知,m為關于N的函數(shù),取k=M-1,則式(23)分母表示為關于N的函數(shù)f(x)：

(24)

對式(23),畫出m關于N的變化曲線如圖4所示,實線,點線,虛線分別為M=64,M=32,M=16時m隨N的變化曲線,可以看到曲線m先增大后減小,又考慮式(18)可知基于卷積的方案運算量與幀的長度無關,所以取在極值點處的幀大小N時,基于重疊相加DTF的優(yōu)化方案的運算量最小。

圖4 m隨N的變化曲線

這個極值點就是式(24)的解,設為x*。所以M一定,基于重疊相加DFT優(yōu)化算法取一幀大小N=[x*]時,運算量最小。

對于已知M和N的情況,顯然當M和N滿足式(25)時基于重疊相加DFT優(yōu)化算法的運算量比較小。

(25)

4 數(shù)值與實驗分析

4.1 MATLAB仿真分析

為了便于對效果的觀察,先假設頻率計權為一個具有低通效果的頻率計權曲線,截止頻率為fs=7 680 Hz。在MATLAB上對基于重疊相加DFT的頻率計權優(yōu)化方案進行仿真,為滿足高精度要求我們選取頻率計權濾波器系數(shù)長度M=2 049,依式(24)得此時最佳的幀大小為：N=28 229,考慮人耳能聽到的聲壓頻率范圍為：20～20 000 Hz,結合采樣定理,選擇采樣頻率f=51 200 Hz,讓三列頻率分別為f1=1 000 Hz,f2=11 000 Hz和f2=17 000 Hz,幅值為1的正弦信號疊加即：x(n)=sin(2πnf1/f)+sin(2πnf2/f)+sin(2πnf3/f),如圖5(a)所示;進行頻率計權得到結果如圖5(b)所示,可以看到高頻成分由于頻率計權增益為0,所以基本得到了抑制,只留下低頻段頻率f1=1 000 Hz幅值為1的正弦信號。

現(xiàn)在對于A計權,讓一列噪聲信號分別通過MATLAB的simulink模塊和使用重疊相加法DFT的方法進行頻率計權,選擇信號采樣頻率f=51 200 Hz,濾波器系數(shù)長度M=1 025,依式(24)得此時最佳的幀大小為：N=13 351,對一列噪聲信號采用simulink仿真和基于重疊相加DFT得結果如圖6(a)所示,實線的是采用simulink模塊頻率計權的結果,點線是使用基于重疊相加DFT法在MATLAB上仿真的結果,可以看到兩條曲線基本重合,也就是兩種方法的計權結果是一致的。為便于觀察將其錯開一定的相位如圖6(b)所示,可以看到兩者的結果是一致的。

圖5 使用DFT計算低通計權結果

圖6 MATLAB simulink 仿真結果

4.2 嵌入式平臺實驗分析

在嵌入式平臺上進行試驗分析[12],驗證算法的準確性,可行性。圖7是實驗平臺實物圖。本文采用DSP+ARM+FPGA的硬件架構,DSP與ARM集成在同一塊雙核處理器上。選擇信號采樣頻率f=51 200 Hz,FPGA負責與慢速外設的交互,聲壓信號AD轉換后,經(jīng)FPGA采集通過SPI傳給DSP，DSP使用EDMA采集得一幀數(shù)據(jù)N=2 048,后進入中斷,這里選擇一幀大小為N=2 048是綜合分析了精度要求和更新速度,系統(tǒng)內存等因素后的結果。為滿足精度要求選擇頻率計權濾波器系數(shù)M=2 049。M和N都已知,代入式(21)得m=26.8,顯然滿足式(25),此時使用基于重疊相加DFT的算法需要的運算量只為基于卷積的方法的運算量的 。DSP處理結果通過DSPLINK傳給ARM,進行實時的分析與顯示,并通過網(wǎng)絡接口導出在MATLAB上進行對比分析。

對同一列噪聲信號分別使用MATLAB仿真和嵌入式平臺實驗得到頻率計權結果如圖8(a)所示,實線是MATLAB仿真的結果,點線是在嵌入式實驗平臺上用重疊相加DFT法計算的到的結果,可以看到兩條曲線基本重合。圖8(b)是把兩條曲線錯開一定相位并拉開后的圖形。也可以看到兩者結果是一致的。

圖7 實驗平臺

圖8 嵌入式系統(tǒng)實驗結果

5 結 論

本文系統(tǒng)分析了基于卷積的頻率計權實現(xiàn)方法和實現(xiàn)過程。但隨著濾波精度的提高,濾波器系數(shù)的增長,計權過程運算量也急劇增大。針對這個問題提出了基于重疊相加DFT的頻率計權優(yōu)化方案,對整個實現(xiàn)過程進行分析,對比和實驗驗證有如下結論：

(1)基于重疊相加DFT的實現(xiàn)方法和基于卷積的實現(xiàn)方法的結果是一致的。

(2)濾波器系數(shù)越長,基于重疊相加DFT的實現(xiàn)方法相比基于卷積的實現(xiàn)方法性能越好。

(3)當濾波器系數(shù)長度確定，x*為方程式(24)的解,取一幀大小N=[x*]時,基于重疊相加DFT的實現(xiàn)方法的運算量最小。

(4)當濾波器系數(shù)長度M,一幀長度N都已知時。若不等式(25)成立,選擇基于重疊相加DFT的實現(xiàn)方法的運算量較小。反之基于卷積的實現(xiàn)方法的運算量較小。

[1]Christopher W T,Bom J K,Chiemi T.Frequency-weighting functions for broadband speech as estimated by a correlational method[J].Acoust.Soc.Am,1998,104(3):1580-1585.

[2]Chong K S，Gwee B H.A 16-channel low-power nonuniform spaced filter bank core for digital hearing aids[J].IEEE Trans.Circuits Syst,2006,53(9):853-857.

[3]楊昌棋,秦樹人,張躍俊.虛擬式噪聲分析儀的數(shù)字計權與開發(fā) [J].重慶大學學報(自然科學版),2001,24(5)：59-66.

YANG Chang-qi,QIN Shu-ren,ZHANG Yue-jun.Digit weight and development of a virtual noise analyzer[J].Journal of Chongqing University(Natural Science Edition),2001,24(5): 59-66.

[4]金暉,何潔.頻率計權的全數(shù)字實現(xiàn) [J].儀器儀表學報,2006,27(6)：1495-1499.

JIN Hui，HE Jie.Digital design method of the frequency weighting[J].Chinese Journal of Scientific Instrument,2006,27(6): 1495-1499.

[5]Zuo L,Nayfeh S A.Low order continuous-time filters for approximation of the ISO 2631-1 human vibration sensitivity weightings[J].Journal of Sound and Vibration,2003,265: 459-465.

[6]Andrew N R,Neil J M.Design of digital filters for frequency weightings required for risk and assessments of workers exposed to vibration [J].Industrial Health，2007,45：512-519.

[7]Emmanuel C L,Barrie W J.Digital signal processing: A practical approach[M].Prentice Hall,2002.

[8]Agarwal R C,Cooley J W.New algorithms for digital convolution[J].IEEE Trans.On ASSP,1977,25(8):281-294.

[9]Steohen A,Martucci.Symmetric convolution and the discrete sine and cosine transforms [J].IEEE Trans.On Signal Processing，1994,42(5)：1038-1051.

[10]胡廣書.數(shù)字信號處理理論,算法與實現(xiàn)(第二版)[M].北京：清華大學出版社，2003.

[11]Oppenheim A V,Schafer R W.Discrete-time signal processing[M].Prentice Hall,1999.

[12]Gajski D D,Vahid F.Specification and design of embedded hardware-software Systems[J].IEEE Des.Test Comput,1995,12(1):53-67.