基于嵌入式系统的频率计权实现方法分析与对比

罗自荣，陈章位，蔡德威，樊亚容

(浙江大学 流体传动与机电系统国家重点实验室,杭州 310013)

近年来人们对噪声污染的关注不断加大,声学分析技术在噪声评估,故障诊断,污染防治,消噪减噪,疾病治疗等方面的应用也不断加深[1-2]。人耳对声音的感受不仅跟声压大小有关,而且跟声压的频率有关。在声学分析中使用频率计权网络,以达到模拟人耳听觉的效果。

模拟电路时代,频率计权的实现主要是通过利用电容,电阻,运放电路等组成计权网络实现频率计权,但是这种方法由于各种电子器件的精度较难保证,而且设计比较复杂,不易变更和调试,同时可靠性和稳定性都比较低。近年来全数字式的频率计权方法得到广泛的应用[3],通过设计一个频率计权滤波器,然后把采样后的声压信号与滤波器做卷积运算,得到计权结果,但随着滤波精度的提高,滤波器系数的增长,计权运算需要的运算量也急剧增大。本文提出了基于重叠相加DFT的优化方案,并对实现过程进行理论分析。对基于卷积和基于重叠相加DFT的方案的运算量级进行理论推导,计算和讨论。理论推导和实验结果表明：两种方案的结果是一致的,随着滤波器系数增大后者的性能快速提升,并得出在滤波器系数已知,和滤波器系数,一帧大小都已知的情况下的方案和参数选择准则。

1 基于卷积的全数字式频率计权实现

1.1 频率计权声压

在声学测量仪器中,依据不同需求使用不同的计权网络,现在常用的频率计权网络有A,B,C 计权。

图1 A,B,C计权网络衰减曲线

图1所示的为A,B和C计权网络衰减曲线图。可以看到人耳对高频段和低频段的声压的反应比较迟钝,对1 000 Hz左右的声压反应比较灵敏。式(1)和(2)分别是A,C计权频率与对应的衰减量的函数关系式。其中：

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 基于卷积的频率计权实现

声压频率计权可通过让声压信号经过对应频率计权滤波器实现。使用双线性变换法,结合(1)式得A,C计权所对应的频率计权滤波器的传递函数[4]

(2)

对频率计权滤波器传递函数从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。

得到滤波器传递函数Z变换,容易得滤波器的系数,假设取滤波器的系数长度为M,调用卷积计算式(3)即可以求得对应频率计权的结果[6]。

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

(4)

如此一帧一帧处理,便可得频率计权后的信号。但是随着分析精度要求的提高,逐渐要求滤波器系数越来越长,而计权滤波的运算量跟滤波器长度有关。容易得到对于一列长为L的数据,通过一个系数长度为M的滤波器总共需要L×M次乘除法和L×(M-1)次加减法运算。所以当M增大时,处理需要的运算量也急剧的增大。

2 基于重叠相加DFT的频率计权实现

2.1 基于重叠相加DFT的卷积实现

随着计权滤波器系数长度的增大,运算量急剧上升,考虑用循环卷积计算线性卷积,而依据卷积定理可通过DFT/IDFT计算循环卷积[7-9]。设滤波器系数h(n)长度为M,一帧长度为N的数据x(n)(x(n)是N个时域声压信号以N为周期的周期延拓,是确定信号,可以使用DFT进行分析),yc(n)和y(n)分别表示h(n)与x(n)的L点循环卷积和线性卷积,且满足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)

联立式(5),式(6)和式(7)得：

(8)

所以h(n)和x(n)的L点循环卷积为线性卷积y(n)以L为周期的周期延拓序列的主值序列。由定义可知y(n)序列的长度为N+M-1。所以y(n)以L为周期的周期延拓无混叠的条件是:

L≥N+M-1

(9)

即当循环卷积长度L,一帧长度N和滤波器系数长度M满足式(5)时有：yc(n)=y(n)。此时依卷积定理通过DFT/IDFT计算得h(n)和x(n)的L点循环卷积,进而得到线性卷积值。具体计算过程如图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)经过频率计权后的结果y(n)为各帧卷积重叠相加。当L=N+M-1时,即为第(i-1)帧的末尾M-1个值与第i帧的前面(M-1)个值相加。具体过程如图3所示。

图3 重叠相加法计算卷积步骤

2.2 频率计权中DFT的实现

采用FFT计算DFT很大得减少了DFT的运算量,但是FFT运算是针对复数的,而要处理的声压信号都是实数,实数可以看为虚部为0的复数,但很明显这有很多的冗余运算,所以考虑把两个相邻的实数看为一个复数[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)联立有：

(14)

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

2.3 频率计权中IDFT的实现

按图2计算得Yc(k),需要进行L点IDFT才能求得最终的频率计权结果。类似DFT的运算,是否有一种简便的方法减少IDFT需要计算的点数。考虑频率计权的声压信号是实数,而设计的滤波器的系数也都是实数,实数与实数的卷积的结果肯定也是实数。也就是说在对Yc(k)进行L点IDFT后的结果是为L个实数。即Yc(k)为yc(n)的L点实序列DFT,设Y1(k)和Y2(k)分别为yc(n)的奇序列和偶序列DFT。则有[11]：

(15)

(16)

另一方面，yc(n)的派生复数序列,即把相邻的两个实数,一个作为复数的实部,一个作为复数的虚部，g(n)=y1(n)+jy2(n),即G(k)=Y1(k)+jY2(k),则联立式(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 算法复杂度分析与对比

设嵌入式系统做一个乘除法运算的时间为Tm,一个加减法运算需要的时间为Ta。滤波器的系数的长度为M,考虑嵌入式系统做卷积实现的需要和避免频繁的数据移动,一般一帧的长度N-1应大于或等于M,则对于直接采用卷积定义的频率计权实现方法计算一列长为Le(Le→∞)的序列x(n)需要的运算量为

(18)

对于基于重叠相加DTF的优化方案,设一帧的长度为N,Le=nN,取循环卷积长度L=N+M-1,则需要的运算量为

(19)

在DSP中采用了硬件乘法器的结构,乘法运算的运算量大大减小,跟加法指令的运算周期接近,这里取Ta=Tm=T,所以：

(20)

对式(20)分子分母同时除以M得：

对一帧大小已知的情况,即N为已知量,m为关于M的函数,显然式(21)的分子为关于变量M的单调增函数。对于式(21)的分母,取k=N-1,则分母可表示为以M为自变量的函数f(x)：

对f(x)求导得：

(22)

考虑前面分析N-1≥M,即x≤k,且滤波器系数M大于1,即：x>1,所以对于式(22)有：

所以f(x)为一个单调减函数,对于式(21),分子是M的单调增函数,分母是M的单调减函数,所以M越大,基于重叠相加DFT优化算法的优越性越好。

对于式(20),分子分母同时除以N得:

(23)

对于滤波器设计已经完成,即M已知,m为关于N的函数,取k=M-1,则式(23)分母表示为关于N的函数f(x)：

(24)

对式(23),画出m关于N的变化曲线如图4所示,实线,点线,虚线分别为M=64,M=32,M=16时m随N的变化曲线,可以看到曲线m先增大后减小,又考虑式(18)可知基于卷积的方案运算量与帧的长度无关,所以取在极值点处的帧大小N时,基于重叠相加DTF的优化方案的运算量最小。

图4 m随N的变化曲线

这个极值点就是式(24)的解,设为x*。所以M一定,基于重叠相加DFT优化算法取一帧大小N=[x*]时,运算量最小。

对于已知M和N的情况,显然当M和N满足式(25)时基于重叠相加DFT优化算法的运算量比较小。

(25)

4 数值与实验分析

4.1 MATLAB仿真分析

为了便于对效果的观察,先假设频率计权为一个具有低通效果的频率计权曲线,截止频率为fs=7 680 Hz。在MATLAB上对基于重叠相加DFT的频率计权优化方案进行仿真,为满足高精度要求我们选取频率计权滤波器系数长度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的正弦信号。

现在对于A计权,让一列噪声信号分别通过MATLAB的simulink模块和使用重叠相加法DFT的方法进行频率计权,选择信号采样频率f=51 200 Hz,滤波器系数长度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转换后,经FPGA采集通过SPI传给DSP，DSP使用EDMA采集得一帧数据N=2 048,后进入中断,这里选择一帧大小为N=2 048是综合分析了精度要求和更新速度,系统内存等因素后的结果。为满足精度要求选择频率计权滤波器系数M=2 049。M和N都已知,代入式(21)得m=26.8,显然满足式(25),此时使用基于重叠相加DFT的算法需要的运算量只为基于卷积的方法的运算量的 。DSP处理结果通过DSPLINK传给ARM,进行实时的分析与显示,并通过网络接口导出在MATLAB上进行对比分析。

对同一列噪声信号分别使用MATLAB仿真和嵌入式平台实验得到频率计权结果如图8(a)所示,实线是MATLAB仿真的结果,点线是在嵌入式实验平台上用重叠相加DFT法计算的到的结果,可以看到两条曲线基本重合。图8(b)是把两条曲线错开一定相位并拉开后的图形。也可以看到两者结果是一致的。

图7 实验平台

图8 嵌入式系统实验结果

5 结 论

本文系统分析了基于卷积的频率计权实现方法和实现过程。但随着滤波精度的提高,滤波器系数的增长,计权过程运算量也急剧增大。针对这个问题提出了基于重叠相加DFT的频率计权优化方案,对整个实现过程进行分析,对比和实验验证有如下结论：

(1)基于重叠相加DFT的实现方法和基于卷积的实现方法的结果是一致的。

(2)滤波器系数越长,基于重叠相加DFT的实现方法相比基于卷积的实现方法性能越好。

(3)当滤波器系数长度确定，x*为方程式(24)的解,取一帧大小N=[x*]时,基于重叠相加DFT的实现方法的运算量最小。

(4)当滤波器系数长度M,一帧长度N都已知时。若不等式(25)成立,选择基于重叠相加DFT的实现方法的运算量较小。反之基于卷积的实现方法的运算量较小。

[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]杨昌棋,秦树人,张跃俊.虚拟式噪声分析仪的数字计权与开发 [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]金晖,何洁.频率计权的全数字实现 [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]胡广书.数字信号处理理论,算法与实现(第二版)[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.