A-Weighting(A计权网络)

前言

通常,所有SLM(声级计)设备都采用频率计权滤波器来调整所测得的噪声值, 而现在市面上流通的声级计一般会提供A和C两种计权的支持,但是通常情况下,我们都采用A计权进行测量,因为A计权与人类的主观感知有很好的相关性,针对人耳对低频信号不那么敏感的情况进行了补偿,并且常常用于测量环境噪声,其测量单位为dBA。

A计权频域实现

在国家标准GB/T 3785.1-2010中,给出了A频率计权滤波器的计算公式,见下式1
A ( f ) = 20 l o g 10 [ f 4 2 f 2 ( f 2 + f 1 2 ) ( f 2 + f 2 2 ) 1 / 2 ( f 2 + f 3 2 ) 1 / 2 ( f 2 + f 4 2 ) ] − A 1 000 (公式1) A(f)=20log_{10}[\frac{f_4^2f^2} { (f^2+f_1^2) (f^2+f_2^2)^{1/2} (f^2+f_3^2)^{1/2} (f^2+f_4^2) } ]-A_1000\tag{公式1} A(f)=20log10[(f2+f12)(f2+f22)1/2(f2+f32)1/2(f2+f42)f42f2]A1000(1)
其中, f 1 = 20.60 H z , f 2 = 107.7 H z , f 3 = 737.9 H z , f 4 = 12194 H z , A 1000 f_1=20.60 Hz,f_2=107.7 Hz,f_3=737.9 Hz,f_4=12194 Hz,A_{1000} f1=20.60Hzf2=107.7Hzf3=737.9Hzf4=12194HzA1000是以分贝表示的归一化常数,相应于在1kHz提供0dB频率计权所需的电增益,修约到最近的0.001dB则有 A 1 000 = − 2.000 d B A_1000=-2.000 dB A1000=2.000dB
其MATLAB代码如下所示。

clc;clear all;
f1=20.60;
f2=107.7;
f3=737.9;
f4=12194;
SampleRate=44100;
A= zeros(1,SampleRate);
dBA= zeros(1,SampleRate);
k= zeros(1,SampleRate);
for i=1:SampleRate;
    k(i)=i;
    f=i^2;
    A(i)=(1.2588966*(148693636*f)/((f+424.3600)*((f+1.1599e+04)*(f+5.4450e+05))^(1/2)*(f+148693636)));
    dBA(i)=20*log10(1.2588966*(148693636*f)/((f+424.3600)*((f+1.1599e+04)*(f+5.4450e+05))^(1/2)*(f+148693636)));
end
A=A(:,11:SampleRate);
dBA=dBA(:,11:SampleRate);
k=k(:,11:SampleRate);
figure();
plot(k,A);
xlabel('Frequency Hz');
ylabel('Gain Amplitude');
figure();
semilogx(k,dBA,'r');
grid on %标注格栅
xlabel('Frequency (Hz)');
ylabel('Gain dB');

由上述代码可得下面的图1跟图2。
图1是在10 Hz - 20 kHz频率范围内的A计权网络幅度谱图。
在这里插入图片描述
图2是图1的对数形式,单位为dB(A)。
在这里插入图片描述

A计权时域实现

在计算A计权瞬时声压级时,较之频域滤波器,时域滤波器不需要经过FFT和IFFT的计算。根据式1,将其dB增益化为幅度增益,即转化为未取对数的公式,并且将自然频率转化为模拟角频率,可得下式2
A ( Ω ) = 1 0 − A 1000 20 Ω 4 2 Ω 4 ( Ω 2 + Ω 1 2 ) ( Ω 2 + Ω 2 2 ) 1 2 ( Ω 2 + Ω 3 2 ) 1 2 ( Ω 2 + Ω 4 2 ) (公式2) A(Ω)=10^{\frac{-A_{1000}}{20}} \frac{Ω_4^2Ω^4} {(Ω^2+Ω_1^2 )(Ω^2+Ω_2^2 )^\frac12 (Ω^2+Ω_3^2 )^\frac12(Ω^2+Ω_4^2 )}\tag{公式2} A(Ω)=1020A1000(Ω2+Ω12)(Ω2+Ω22)21(Ω2+Ω32)21(Ω2+Ω42)Ω42Ω4(2)
为了设计数字滤波器,我们可以先通过上式推出S复数域的传递函数,然后利用双线性变换法设计对应的IIR数字滤波器。我们首先来观察一个一阶系统,如下式3所示。
H ( S ) = 1 S + Ω 1 (公式3) H(S)=\frac1{S+Ω_1}\tag{公式3} H(S)=S+Ω11(3)
由信号与系统的知识,令公式3中的S=JΩ,可以得到系统的频域表示,并求其系统的幅频响应,如下式4所示。
∣ H ( S ) ∣ = 1 Ω 2 + Ω 1 2 (公式4) |H(S)|=\sqrt{\frac1{Ω^2+Ω_1^2}}\tag{公式4} H(S)=Ω2+Ω121 (4)
通过分析,我们发现式4很像式2中分母的一部分,而一个系统可以看成多个子系统级联起来的,级联的幅频响应是乘积关系,因此我们很容易联想到公式2是由6个一阶系统、1个四阶微分系统和常数增益环节构成,因此级联后的A计权系统的传递函数如下式5所示。
A ( S ) = 1 0 − A 1000 20 Ω 4 2 S 4 ( S 2 + Ω 1 2 ) 2 ( S 2 + Ω 2 2 ) ( S 2 + Ω 3 2 ) ( S 2 + Ω 4 2 ) 2 (公式5) A(S)=10^ \frac{-A_{1000}}{20} \frac{Ω_4^2S^4} {(S^2+Ω_1^2 )^2 (S^2+Ω_2^2 )(S^2+Ω_3^2 ) (S^2+Ω_4^2 )^2 }\tag{公式5} A(S)=1020A1000(S2+Ω12)2(S2+Ω22)(S2+Ω32)(S2+Ω42)2Ω42S4(5)
接下来利用MATLAB里面的双线性变换函数,计算滤波器系数,代码如下所示。

clc;clear all;
Fs=44100;
f1 = 20.598997; 
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
A1000 =1.9997;
pi = 3.14159265358979;
NUMs = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
DENs = conv([1 4*pi*f4 (2*pi*f4)^2],[1 4*pi*f1 (2*pi*f1)^2]); 
DENs = conv(conv(DENs,[1 2*pi*f3]),[1 2*pi*f2]);
[B,A] = bilinear(NUMs,DENs,Fs);
fvtool(B,A)                     % Visualize the filter

通过上述的MATLAB代码,即可得到该数字滤波器的分子分母系数数组分别如下所示( f s = 44100 H z f_s=44100Hz fs=44100Hz)。

%****************** the runing result **************************
% B = 0.2557   -0.5115   -0.2557    1.0230   -0.2557   -0.5115    0.2557
% A = 1.0000   -4.0196    6.1894   -4.4532    1.4208   -0.1418    0.0044
%******************* the runing result**************************

其滤波器的幅频响应如下图3所示。
在这里插入图片描述
我们采用IIR直接Ⅱ型来实现上述滤波器,其公式如下式6所示。
y ( n ) = ∑ i = 0 N b i x ( n − i ) − ∑ j = 1 N a j y ( n − j ) (公式6) y(n)=\sum_{i=0}^N{b_ix(n-i)}-\sum_{j=1}^N{a_jy(n-j)}\tag{公式6} y(n)=i=0Nbix(ni)j=1Najy(nj)(6)
其中,x为采样点输入数组,y为经过A计权采样点的输出数组,n为数组内采样点编号,且1≤n≤BS(Block Size);a为滤波器分母系数数组,b为滤波器分子系数数组, N为数组下标,且 0≤N≤6;假设所有初试条件为0。式中可以看出, FIR滤波器的输出不仅与现在的输入有关,还与过去的输入有关。
最后,利用代码实现式6,Java代码如下所示。

public static double[] IIRFilter(double[] signal, double[] a,double[] b) {
    	double[]in = new double[b.length];
    	double[]out = new double[a.length-1];
    	double[]outData = new double[signal.length];
        for (int i = 0; i < signal.length; i++) {
            System.arraycopy(in, 0, in, 1, in.length - 1);
            in[0] = signal[i];
            //calculate y based on a and b coefficients
            //and in and out.
            float y = 0;
            for(int j = 0 ; j < b.length ; j++){
                y += b[j] * in[j];
            }
            for(int j = 0;j < a.length-1;j++){
                y -= a[j+1] * out[j];
            }
            //shift the out array
            System.arraycopy(out, 0, out, 1, out.length - 1);
            out[0] = y;
            outData[i] = y;
        }
        return outData;
    }

总结

前面我们讲述了A-Weighting(A计权网络)的频域和时域的原理思路及其实现方式,并贴出了相应的实现代码,希望读者们能各取所需,如有问题,欢迎读者们批评斧正 。

感谢

感谢默默为我们提供优质博客讲解的作者们,非常感谢!
下面附上参考资料的链接:
A计权时域实现原理
IIR滤波器Java实现
如有任何侵权的地方,希望作者联系我,我将移除相关内容,再次感谢!

  • 46
    点赞
  • 105
    收藏
    觉得还不错? 一键收藏
  • 40
    评论
评论 40
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值