src/utils/loudness/loudness.sci

   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
  14 function 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
  26 endfunction
  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
  40 function [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
  47 if %nargin<1 then
  48    fs=44100;
  49 end
  50
  51 // Specify the 80 dB Equal Loudness curve
  52 EL80=[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
  61 for 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
  72 end
  73
  74 // convert frequency and amplitude of the equal loudness curve into format suitable for yulewalk
  75 f=EL80(:,1)./(fs/2);
  76 m=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
  82 hz=iir(2,'hp','butt',[150/fs,0],[1e-3 1e-3]);
  83 b2=numer(hz); // b2=b2($:-1:1);
  84 a2=denom(hz); // a2=a2($:-1:1);
  85
  86 endfunction
  87
  88 // *****************************************************************************
  89 // Generate Filter Taps
  90 //
  91 // Generate the filter taps for each of the desired frequencies.
  92
  93 format('v', 16);
  94
  95 freqs = [  8000 11025 12000 16000 18900 22050 24000 ..
  96           28000 32000 36000 37800 44100 48000 ];
  97
  98 for 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
 112 end
 113
 114
 115 quit