а-взвешивание в качестве цифрового фильтра

Question

Я пытаюсь создать цифровой фильтр, который создает близкие результаты в спецификацию а-весовой: reference

Я уже создал простой низко- и верхние частоты с 2 ^ем и 6 ^th заказать фильтры IIR. Я также создал каскадную последовательность второго порядка из трех биквадров. Эти фильтры, похоже, работают достаточно хорошо. Я генерирую коэффициенты с помощью функций октавы/матлаба butter и tf2sos. Фильтр реализован на C++ и применяется к моему аудиовходу. Обработанное аудио преобразуется в частотную область и отображается на моем экране.

Теперь мой вопрос, как я могу реализовать фильтр a-weighting. Я уже попробовал код 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 

, но эти коэффициенты также не работают. Может быть, потому, что коэффициенты ошибочны в первую очередь?

Может кто-нибудь указать мне, как правильно реализовать цифровой весовой фильтр?

EDIT: Благодаря ответ SleuthEye в поле ниже, я мог бы реализовать а-взвешивающий фильтр с правильными коэффициентами из октавы.

Я нашел множество различных реализаций фильтра C++ в Интернете и в приложениях-генераторах win32-фильтров, но они казались все непоследовательными.

Наконец, я нашел два кодовых фрагментов, которые в основном работают таким же образом:

• 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/ (а и б поменяны местами здесь)

Так что я в основном скопировал код оттуда и закончил с этим: (его транспонированная прямая реализация 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; 
}; 

Initialization:

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's bilinear (известный как incompatibilities with Matlab), для которого третьим аргументом является период выборки (вместо частоты выборки для 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'); 

Это должно дать вам график ниже (синяя кривая является эталонной, а красная - получена из ваших коэффициентов):