a-weighting as digital filter

I'm trying to create a digital filter that creates close results to the a-weighting specification: reference

I have already created simple low- and high-pass with 2^nd and 6^th order IIR filters. I also created a cascaded second order sequence of 3 biquads. These filters seem to work quite well. I'm generating the coefficients with the octave/matlab butter and tf2sos functions. The filter is implemented in C++ and applied to my audio input. The processed audio gets converted to the frequency domain and is displayed on my screen.

Now is my question how I can implement the a-weighting filter. I already tried the matlab code I found on the internet:

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

These coefficients look strange, and they also do not work with my 6^th order IIR. So I tried to create cascaded biquads:

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

but these coefficients do not work either. Could it be because the coefficients are wrong in the first place?

Can someone point me out, how to implement a digital a-weighting filter correctly?

EDIT: Thanks to SleuthEye's answer below, I could implement the a-weighting filter with the correct coefficients from octave.

I found a lot of different filter c++ implementations on the web and in win32-filter-generator applications, but they seemed all inconsistent.

Finally I found two code-snippets that essentially work in the same way:

https://github.com/berndporr/iir1/blob/master/iir/State.h (has also a youtube channel: https://www.youtube.com/user/DSPcourse)
http://www.earlevel.com/main/2012/11/26/biquad-c-source-code/ (a and b are swapped here)

So I basically copied the code from there and ended up with this: (its the transposed direct form II implementation)

Definition:

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

Invocation:

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

Solution

To get the results you quote it would seem that you are actually using Octave's 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).

To work around this problem, you can get the sampling period as 1/Fs and correspondingly compute your coefficients with:

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

You should then be able to confirm the filter response with:

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

This should give you the graph below (blue curve being the reference, and the red one obtained from your coefficients):