// 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