soft/test/audio/std_utils/stdweight.m

87 lines
3.3 KiB
Matlab

function y = stdweight(x, fs)
%% y = stdweight(x, fs)
%
% Input
% x - input signal
% fs - sample rate
%
% Output
% y - filtered signal
%
%%
% Copyright (c) 2017, Intel Corporation
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
% * Neither the name of the Intel Corporation nor the
% names of its contributors may be used to endorse or promote products
% derived from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
%
%% Frequency grid for FIR design
f = linspace(0, fs/2, 6000);
%% From https://en.wikipedia.org/wiki/ITU-R_468_noise_weighting
h1 = -4.737338981378384e-24*f.^6 +2.043828333606125e-15*f.^4 -1.363894795463638e-7*f.^2 +1;
h2 = 1.306612257412824e-19*f.^5 -2.118150887518656e-11*f.^3 +5.559488023498642e-4*f;
RITU = 1.246332637532143e-4*f./sqrt(h1.^2+h2.^2);
ITU = 18.2+20*log10(RITU)-5.63; % AES17 CCIR-RMS
fc = [ 31.5 63 100 200 400 800 1000 2000 3150 4000 5000 6300 7100 8000 9000 10000 12500 14000 16000 20000 31500];
ga = [-35.5 -29.5 -25.4 -19.4 -13.4 -7.5 -5.6 0 3.4 4.9 6.1 6.6 6.4 5.8 4.5 2.5 -5.6 -10.9 -17.3 -27.8 -48.3];
tu = [ 2.0 1.4 1.0 0.85 0.70 0.55 0.50 0.50 0.50 0.50 0.50 0.01 0.20 0.40 0.60 0.80 1.2 1.4 1.6 2.0 2.8];
tl = [ -2.0 -1.4 -1.0 -0.85 -0.70 -0.55 -0.50 -0.50 -0.50 -0.50 -0.50 -0.01 -0.20 -0.40 -0.60 -0.80 -1.2 -1.4 -1.6 -2.0 NaN];
%% Fine tune the 6300Hz, 6.6 dB point
if fs > 2*6300
idx = find(f > 6300, 1, 'first');
ITU = ITU + 6.6-ITU(idx);
else
idx = find(f > 3150, 1, 'first');
ITU = ITU + 3.4-ITU(idx);
end
m_fir = 10.^(ITU/20);
f_fir = 2*f/fs;
n_fir = 4000; % Sufficient up to 192 kHz
if (fs < 96001)
n_fir = 2000;
end
if (fs < 48001)
n_fir = 1000;
end
if (fs < 16001)
n_fir = 400;
end
bz = fir2(n_fir, f_fir, m_fir);
y = filter(bz, 1, x);
if 0
h = freqz(bz, 1, f, fs);
my_w_db = 20*log10(abs(h));
figure;
semilogx(f,ITU,fc,ga,'x',fc,ga+tu,fc,ga+tl, f, my_w_db );
grid on; xlabel('Hz'); ylabel('dB');
end
end