Port David Robinson's equalloudfilt.m MATLAB script to scilab so that the ReplayGain...

Signed-off-by: Earl Chew <earl_chew@yahoo.com>

diff --git a/src/utils/loudness/loudness.sci b/src/utils/loudness/loudness.sci
new file mode 100644 (file)
index 0000000..a476a17
--- /dev/null
+++ b/src/utils/loudness/loudness.sci
@@ -0,0 +1,115 @@
+// 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