function [b, a] = aweighting(fs, type)
%AWEIGHTING Digital A-weighting in time domain
% [b, a] = AWEIGHTING(fs) returns the discrete time realization of the
% A-weighting filter.
% [b, a] = AWEIGHTING(fs, type) returns the general weighting filter, where
% type can be 'A', 'B' or 'C'.

% Copyright (c) 2013-2016 Peter Fiala

if nargin < 2
    type = 'A';
end

% equivalent continuous time poles and zeros from Wikipedia
switch upper(type)
    case 'A'
        % zeros, poles, amplification
        z = [0 0 0 0];
        p = [-129.4 -129.4 -676.7 -4636 -76655 -76655];
        k = 7.3905e9;
    case 'B'
        % zeros, poles, amplification
        z = [0 0 0];
        p = [-129.4 -129.4 -995.9 -76655 -76655];
        k = 5.99185e9;
    case 'C'
        % zeros, poles, amplification
        z = [0 0];
        p = [-129.4 -129.4 -76655 -76655];
        k = 5.91797e9;
    otherwise
        error('aweighting:invalid_argument',...
            'Argument ''type'' must be ''A'', ''B'' or ''C''');
end

% continuous system
filterc = zpk(z, p, k);
% convert to discrete
filterd = c2d(filterc, 1/fs);
% get transfer function equivalent
[num, den] = tfdata(filterd);
b = num{1};
a = den{1};

end % of function aweighting
