mirror of
https://github.com/claunia/flac.git
synced 2025-12-16 18:54:26 +00:00
116 lines
3.1 KiB
Scilab
116 lines
3.1 KiB
Scilab
// Equal Loudness Filter
|
|
//
|
|
// Adapted from original MATLAB code written by David Robinson
|
|
//
|
|
// http://replaygain.hydrogenaudio.org/proposal/equal_loudness.html
|
|
// http://replaygain.hydrogenaudio.org/proposal/mfiles/equalloudfilt.m
|
|
|
|
// *****************************************************************************
|
|
// Print Filter Coefficients
|
|
//
|
|
// This function takes a vector of filter tap settings, and prints
|
|
// each tap setting from least significant to most significant.
|
|
|
|
function c=printcoeff(p)
|
|
|
|
c=coeff(p);
|
|
c=c($:-1:1);
|
|
|
|
for ix = 1:1:length(c)
|
|
if ix > 1
|
|
printf(" ")
|
|
end
|
|
printf("%.14f", c(ix));
|
|
end
|
|
|
|
endfunction
|
|
|
|
// *****************************************************************************
|
|
// Equal Loudness Filter
|
|
//
|
|
// This function is adapted from David Robison's original MATLAB code.
|
|
// Apart from changes to port it to scilab, the other change is to
|
|
// use a single specification of the frequency points in the 80dB Equal
|
|
// Loudness curve.
|
|
//
|
|
// The original code had different curves for different sampling
|
|
// frequencies. This code dynamically computes the current data
|
|
// points to use as determined by the Nyquist frequency.
|
|
|
|
function [a1,b1,a2,b2]=equalloudfilt(fs);
|
|
// Design a filter to match equal loudness curves
|
|
// 9/7/2001
|
|
|
|
[%nargout,%nargin]=argn(0);
|
|
|
|
// If the user hasn't specified a sampling frequency, use the CD default
|
|
if %nargin<1 then
|
|
fs=44100;
|
|
end
|
|
|
|
// Specify the 80 dB Equal Loudness curve
|
|
EL80=[0,120;20,113;30,103;40,97;50,93;60,91;70,89;80,87;90,86; ..
|
|
..
|
|
100,85;200,78;300,76;400,76;500,76;600,76;700,77;800,78;900,79.5; ..
|
|
..
|
|
1000,80;1500,79;2000,77;2500,74;3000,71.5;3700,70;4000,70.5; ..
|
|
5000,74;6000,79;7000,84;8000,86;9000,86; ..
|
|
..
|
|
10000,85;12000,95;15000,110;20000,125;24000,140];
|
|
|
|
for ex = 1:1:length(EL80(:,1))
|
|
if EL80(ex,1) > fs/2
|
|
EL80 = [ EL80(1:ex-1,:); fs/2, EL80(ex-1,2) ];
|
|
break
|
|
elseif EL80(ex,1) == fs/2
|
|
EL80 = EL80(1:ex,:);
|
|
break
|
|
end
|
|
if ex == length(EL80(:,1))
|
|
EL80 = [ EL80(1:$, :); fs/2, EL80($,2) ];
|
|
end
|
|
end
|
|
|
|
// convert frequency and amplitude of the equal loudness curve into format suitable for yulewalk
|
|
f=EL80(:,1)./(fs/2);
|
|
m=10.^((70-EL80(:,2))/20);
|
|
|
|
// Use a MATLAB utility to design a best bit IIR filter
|
|
[b1,a1]=yulewalk(10,f,m);
|
|
|
|
// Add a 2nd order high pass filter at 150Hz to finish the job
|
|
hz=iir(2,'hp','butt',[150/fs,0],[1e-3 1e-3]);
|
|
b2=numer(hz); // b2=b2($:-1:1);
|
|
a2=denom(hz); // a2=a2($:-1:1);
|
|
|
|
endfunction
|
|
|
|
// *****************************************************************************
|
|
// Generate Filter Taps
|
|
//
|
|
// Generate the filter taps for each of the desired frequencies.
|
|
|
|
format('v', 16);
|
|
|
|
freqs = [ 8000 11025 12000 16000 18900 22050 24000 ..
|
|
28000 32000 36000 37800 44100 48000 ];
|
|
|
|
for fx = 1:1:length(freqs)
|
|
|
|
printf("\n%d\n", freqs(fx));
|
|
|
|
[a1,b1,a2,b2] = equalloudfilt(freqs(fx));
|
|
|
|
printf("{ "); bb=printcoeff(b1); printf(" }\n");
|
|
printf("{ "); aa=printcoeff(a1); printf(" }\n");
|
|
|
|
printf("{ "); printcoeff(b2); printf(" }\n");
|
|
printf("{ "); printcoeff(a2); printf(" }\n");
|
|
|
|
// freqz_fwd(bb,aa,1024,freqs(fx));
|
|
|
|
end
|
|
|
|
|
|
quit
|