--- /dev/null
+// 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