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

(浙江大學(xué) 流體傳動與機(jī)電系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室,杭州 310013)

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

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

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

1.1 頻率計(jì)權(quán)聲壓

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

圖1 A,B,C計(jì)權(quán)網(wǎng)絡(luò)衰減曲線

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

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 基于卷積的頻率計(jì)權(quán)實(shí)現(xiàn)

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

(2)

對頻率計(jì)權(quán)濾波器傳遞函數(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è)取濾波器的系數(shù)長度為M,調(diào)用卷積計(jì)算式(3)即可以求得對應(yīng)頻率計(jì)權(quán)的結(jié)果[6]。

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

(4)

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

2 基于重疊相加DFT的頻率計(jì)權(quán)實(shí)現(xiàn)

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

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

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

(5)

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

(6)

又依卷積的定義有：

(7)

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

(8)

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

L≥N+M-1

(9)

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

圖2 卷積計(jì)算流程圖

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

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

(10)

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

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

(11)

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

圖3 重疊相加法計(jì)算卷積步驟

2.2 頻率計(jì)權(quán)中DFT的實(shí)現(xiàn)

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

(12)

又結(jié)合DFT的定義和性質(zhì)有：

(13)

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

(14)

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

2.3 頻率計(jì)權(quán)中IDFT的實(shí)現(xiàn)

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

(15)

(16)

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

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

(17)

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

3 算法復(fù)雜度分析與對比

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

(18)

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

(19)

在DSP中采用了硬件乘法器的結(jié)構(gòu),乘法運(yùn)算的運(yùn)算量大大減小,跟加法指令的運(yùn)算周期接近,這里取Ta=Tm=T,所以：

(20)

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

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

對f(x)求導(dǎo)得：

(22)

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

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

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

(23)

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

(24)

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

圖4 m隨N的變化曲線

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

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

(25)

4 數(shù)值與實(shí)驗(yàn)分析

4.1 MATLAB仿真分析

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

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

圖5 使用DFT計(jì)算低通計(jì)權(quán)結(jié)果

圖6 MATLAB simulink 仿真結(jié)果

4.2 嵌入式平臺實(shí)驗(yàn)分析

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

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

圖7 實(shí)驗(yàn)平臺

圖8 嵌入式系統(tǒng)實(shí)驗(yàn)結(jié)果

5 結(jié) 論

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

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

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

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

(4)當(dāng)濾波器系數(shù)長度M,一幀長度N都已知時。若不等式(25)成立,選擇基于重疊相加DFT的實(shí)現(xiàn)方法的運(yùn)算量較小。反之基于卷積的實(shí)現(xiàn)方法的運(yù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ù)字計(jì)權(quán)與開發(fā) [J].重慶大學(xué)學(xué)報(bào)(自然科學(xué)版),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]金暉,何潔.頻率計(jì)權(quán)的全數(shù)字實(shí)現(xiàn) [J].儀器儀表學(xué)報(bào),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ù)字信號處理理論,算法與實(shí)現(xiàn)(第二版)[M].北京：清華大學(xué)出版社，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.