A-Weighting（A计权网络）

前言

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

A计权频域实现

在国家标准GB/T 3785.1-2010中，给出了A频率计权滤波器的计算公式，见下式1
$$A(f)=20lo{g}_{10}[\frac{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\text{\hspace{1em}(公式1)}$$
其中，$f_1=20.60Hz，f_2=107.7Hz，f_3=737.9Hz，f_4=12194Hz，A_{1000}$是以分贝表示的归一化常数，相应于在1kHz提供0dB频率计权所需的电增益，修约到最近的0.001dB则有$A_{1}000=-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
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

由上述代码可得下面的图1跟图2。
图1是在10 Hz - 20 kHz频率范围内的A计权网络幅度谱图。

图2是图1的对数形式，单位为dB（A）。

A计权时域实现

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

通过上述的MATLAB代码，即可得到该数字滤波器的分子分母系数数组分别如下所示($f_s=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**************************
1
2
3
4

其滤波器的幅频响应如下图3所示。

我们采用IIR直接Ⅱ型来实现上述滤波器，其公式如下式6所示。
$$y(n)=\sum _{i=0}^N{b}_{i}x(n-i)-\sum _{j=1}^N{a}_{j}y(n-j)\text{\hspace{1em}(公式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;
    }
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

总结

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

感谢

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