xref: /aosp_15_r20/external/flac/src/utils/loudness/loudness.sci (revision 600f14f40d737144c998e2ec7a483122d3776fbc)
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