a加权作为数字滤波器

Question

我正在尝试创建一个数字滤波器，以创建接近 a 加权规范的结果：reference

我已经用 2^nd 和 6^th 阶 IIR 滤波器创建了简单的低通和高通。我还创建了 3 个双二阶的级联二阶序列。这些过滤器似乎工作得很好。我正在使用 octave/matlab butter 和 tf2sos 函数生成系数。该过滤器是用 C++ 实现的，并应用于我的音频输入。处理后的音频转换为频域并显示在我的屏幕上。

现在我的问题是如何实现 a 加权滤波器。我已经尝试了我在互联网上找到的 matlab 代码：

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);

printf('%0.16f\n',A);
1.0000000000000000
5.9999984423677502
14.9999922118394622
19.9999844236803490
14.9999844236817736
5.9999922118415991
0.9999984423684625

printf('%0.16f\n',B);
0.0000000000000000
-0.0000000000000000
-0.0000000000000000
0.0000000000000000
-0.0000000000000000
-0.0000000000000000
0.0000000000000000

这些系数看起来很奇怪，而且它们也不适用于我的第 6 阶 IIR。所以我尝试创建级联双二阶函数：

sos = tf2sos(B,A);
printf('%0.16f\n',sos);
0.0000000000000000
1.0000000000000000
1.0000000000000000
-0.0000000000000000
-2.0001419687327253
1.9999999999999987
0.0000000000000000
1.0001419788113322
0.9999999999999984
1.0000000000000000
1.0000000000000000
1.0000000000000000
2.0044328398652649
1.9955456521230652
2.0000199503794249
1.0044526299594496
0.9955653535664766
1.0000002046824734

但是这些系数也不起作用。会不会是系数一开始就错了？

谁能指出我如何正确实现数字 a 加权滤波器？

编辑： 感谢下面 SleuthEye 的回答，我可以使用来自八度的正确系数实现 a 加权滤波器。

我在网络上和 win32-filter-generator 应用程序中发现了很多不同的过滤器 c++ 实现，但它们似乎都不一致。

最后我找到了两个基本上以相同方式工作的代码片段：

https://github.com/berndporr/iir1/blob/master/iir/State.h （还有一个 youtube 频道：https://www.youtube.com/user/DSPcourse）
http://www.earlevel.com/main/2012/11/26/biquad-c-source-code/（这里a和b互换）

所以我基本上是从那里复制代码并最终得到这个：（它是转置直接形式 II 实现）

定义：

class Biquad{

public:

    Biquad(double a0, double a1, double a2, double b0, double b1, double b2) {
        this->a0 = a0;
        this->a1 = a1;
        this->a2 = a2;
        this->b0 = b0;
        this->b1 = b1;
        this->b2 = b2;
        z1 = z2 = 0.0;
    }

    double filter(double in) {
        double out = z1 + b0 * in;
        z1 = z2 + b1 * in - a1 * out;
        z2 = b2 * in - a2 * out;
        return out;
    }

private:
    double a0, a1, a2, b0, b1, b2;
    double z1, z2;
};

初始化：

Biquad *a1Filter = new Biquad(1.0000000, -0.1405361, 0.0049376, 0.2557411, -0.5114387, 0.2556976);
Biquad *a2Filter = new Biquad(1.0000000, -1.8849012, 0.8864215, 1.0000000, -2.0001702, 1.0001702);
Biquad *a3Filter = new Biquad(1.0000000, -1.9941389, 0.9941475, 1.0000000, 2.0000000, 1.0000000);

调用：

double outputValue = a3Filter->filter(a2Filter->filter(a1Filter->filter(inputValue)));

Answer 1

要获得您引用的结果，您似乎实际上是在使用 Octave 的 bilinear (known to have some incompatibilities with Matlab) for which the third argument is the sampling period (instead of the sampling frequency for Matlab's bilinear)。

要解决此问题，您可以将采样周期设置为 1/Fs 并相应地计算您的系数：

[B,A] = bilinear(NUMs,DENs,1/Fs);

然后您应该能够通过以下方式确认过滤器响应：

fmin = 10;    % Hz
fmax = 20000; % Hz
omega = logspace(log10(2*pi*fmin), log10(2*pi*fmax));
Ha = freqs(NUMs,DENs,omega);
hold off; semilogx(omega/(2*pi), 20*log10(abs(Ha)), 'b');

N = ceil(Fs/fmin);
[Hd,W] = freqz(B,A,N);
hold on; semilogx(W*Fs/(2*pi), 20*log10(abs(Hd)), 'r');

这应该给你下面的图表（蓝色曲线是参考，红色曲线是从你的系数中获得的）：

a加权作为数字滤波器

a-weighting as digital filter

matlab

signal-processing

digital-filter