1// Equal Loudness Filter 2// 3// Adapted from original MATLAB code written by David Robinson 4// 5// http://replaygain.hydrogenaudio.org/proposal/equal_loudness.html 6// http://replaygain.hydrogenaudio.org/proposal/mfiles/equalloudfilt.m 7 8// ***************************************************************************** 9// Print Filter Coefficients 10// 11// This function takes a vector of filter tap settings, and prints 12// each tap setting from least significant to most significant. 13 14function c=printcoeff(p) 15 16 c=coeff(p); 17 c=c($:-1:1); 18 19 for ix = 1:1:length(c) 20 if ix > 1 21 printf(" ") 22 end 23 printf("%.14f", c(ix)); 24 end 25 26endfunction 27 28// ***************************************************************************** 29// Equal Loudness Filter 30// 31// This function is adapted from David Robison's original MATLAB code. 32// Apart from changes to port it to scilab, the other change is to 33// use a single specification of the frequency points in the 80dB Equal 34// Loudness curve. 35// 36// The original code had different curves for different sampling 37// frequencies. This code dynamically computes the current data 38// points to use as determined by the Nyquist frequency. 39 40function [a1,b1,a2,b2]=equalloudfilt(fs); 41// Design a filter to match equal loudness curves 42// 9/7/2001 43 44[%nargout,%nargin]=argn(0); 45 46// If the user hasn't specified a sampling frequency, use the CD default 47if %nargin<1 then 48 fs=44100; 49end 50 51// Specify the 80 dB Equal Loudness curve 52EL80=[0,120;20,113;30,103;40,97;50,93;60,91;70,89;80,87;90,86; .. 53 .. 54 100,85;200,78;300,76;400,76;500,76;600,76;700,77;800,78;900,79.5; .. 55 .. 56 1000,80;1500,79;2000,77;2500,74;3000,71.5;3700,70;4000,70.5; .. 57 5000,74;6000,79;7000,84;8000,86;9000,86; .. 58 .. 59 10000,85;12000,95;15000,110;20000,125;24000,140]; 60 61for ex = 1:1:length(EL80(:,1)) 62 if EL80(ex,1) > fs/2 63 EL80 = [ EL80(1:ex-1,:); fs/2, EL80(ex-1,2) ]; 64 break 65 elseif EL80(ex,1) == fs/2 66 EL80 = EL80(1:ex,:); 67 break 68 end 69 if ex == length(EL80(:,1)) 70 EL80 = [ EL80(1:$, :); fs/2, EL80($,2) ]; 71 end 72end 73 74// convert frequency and amplitude of the equal loudness curve into format suitable for yulewalk 75f=EL80(:,1)./(fs/2); 76m=10.^((70-EL80(:,2))/20); 77 78// Use a MATLAB utility to design a best bit IIR filter 79[b1,a1]=yulewalk(10,f,m); 80 81// Add a 2nd order high pass filter at 150Hz to finish the job 82hz=iir(2,'hp','butt',[150/fs,0],[1e-3 1e-3]); 83b2=numer(hz); // b2=b2($:-1:1); 84a2=denom(hz); // a2=a2($:-1:1); 85 86endfunction 87 88// ***************************************************************************** 89// Generate Filter Taps 90// 91// Generate the filter taps for each of the desired frequencies. 92 93format('v', 16); 94 95freqs = [ 8000 11025 12000 16000 18900 22050 24000 .. 96 28000 32000 36000 37800 44100 48000 ]; 97 98for fx = 1:1:length(freqs) 99 100 printf("\n%d\n", freqs(fx)); 101 102 [a1,b1,a2,b2] = equalloudfilt(freqs(fx)); 103 104 printf("{ "); bb=printcoeff(b1); printf(" }\n"); 105 printf("{ "); aa=printcoeff(a1); printf(" }\n"); 106 107 printf("{ "); printcoeff(b2); printf(" }\n"); 108 printf("{ "); printcoeff(a2); printf(" }\n"); 109 110// freqz_fwd(bb,aa,1024,freqs(fx)); 111 112end 113 114 115quit 116