
// Critical bandwidth vs. frequency
//
// This is the fundamentl starting point for the critical 
// band rate scale as it corresponds to the data gathered
// from different experiments involving masking and other  
// psychoacoustical properties of hearing
// 
// References
// 
// Zwicker (1961) "Subdivision of the Audible Frequency
//     Range into Critical Bands (Frequenzgruppen)"
//     JASA 33 No 2, p. 248.
// Zwicker, Flottorp, Stevens (1957) "Critical Band Width
//     in Loudness Summation". JASA 29 No 5, pp. 548-557.
//
// ---------------------------------
// Federico Miyara
// 2019-04-10
// 2024-07-22

// 1.
// The following data correspond to the standard table
// by Zwicker (1961). As the author comments, this series
// of adjacent critical bands have been chosen in such a
// fashion as to ensure that some frequencies, such as
// 250 Hz, 1000 Hz and 4000 Hz, which are the centers of
// some of the standard octave and one third octave 
// frequency bands, be also centers of some of the
// critical bands. To that end, some rounding of the
// original data had to be accepted. The rounding process
// is not explained and the raw data are not available.
f_zw =   [  50  150  250  350   450   570  700  840 1000 ..
          1170 1370 1600 1850  2150  2500 2900 3400 4000 ..
          4800 5800 7000 8500 10500 13500];

cbw_zw = [  80  100  100  100   110   120  140  150  160 ..
           190  210  240  280   320   380  450  550  700 ..
           900 1100 1300 1800  2500  3500];

// 2.
// Data from Zwicker, Flottorp, Stevens (1957) "Critical
// Band Width in Loudness Summation". JASA 29, 548-557.
// These data come from two sources in the referenced
// paper: Table 1 and Figure 12. Data from Table 1 have
// not been rounded to get some frequencies from the 
// standard one third octave bands, as has been done
// with the most widely used critical bandewidth table, 
// but have certain irregularities. Data from Figure 12
// are a fitting for a limited frequency range.

// Center-left column of table I (column 2)
f1 = [   20  110  200  295   395   503   625  755  900 ..
       1060 1250 1460 1700  1970  2290  2670 3120 3680 ..
       4360 5200 6200 7500  9300 11700 15000];

// Center-right column of table I (column 3)
f2 = [   65  155  250  345   450   560   690  830  980 ..
       1155 1355 1580 1835  2130  2480  2900 3400 4020 ..
       4780 5700 6850 8400 10500 13300 17300];

// As each series of frequencies have been chosen so that
// they differ in one critical band, critical bands are 
// just the difference between two consecutive frequencies
// (they are also included in the table as  columns 1 and
// 4) 
cbw1 = f1(2:$) - f1(1:$-1);
cbw2 = f2(2:$) - f2(1:$-1);

// Both interleaved series of frequencies are put
// consecutively and sorted in ascending order. The 
// original position of the frequencies is saved in
// vector of indices i. 
[f, i] = gsort([f1, f2], 'g', 'i');
// The first and last ones are removed since no critical
// band corresponds to them 
f = f(2:$-1);
i = i(2:$-1);

// Both series of critical bands are put consecutively.
// However, since the first remaining index is 26 (instead
// of 2), corresponding to the second series, first comes
// cbw2 and then cbw1
cbw = [cbw2, cbw1];
// Critical bands are sorted to match the corresponding
// frequencies
cbw = cbw(i-1);
// NOTE: It would have been possible to simply sort cbw 
//       in ascending order as well: 
//       cbw = gsort([cbw2,cbw1], 'g', 'i')
//       but it is not a robust approach since it relies 
//       on the fact that critical bands increase
//       monotonically.    

// 3.
// Raw, manually parametrized data from figure 12 of
// paper. The relationship between x, y and f, cbw is
// x = 10*log10(f/100)
// y = 10*log10(cbw/10)
x = [    0     1     2     3     4     5     6     7     8 ..
         9    10    11    12    13    14    15    16    17 ..
        18    19    20 20.75];
y = [  9.6   9.6  9.63   9.7  9.75  9.92 10.18 10.52  11.0 ..
     11.58 12.25  13.0 13.88  14.9 15.95  17.0  18.2 19.43 ..
      20.7 22.18 23.65 24.85];

// Retrieval of frequency and critical bandwidth
f0 = 100*10.^(x/10);
cbw0 =  10*10.^(y/10);


// Plot the data
scf(1);
clf(1);
// Plot data from table I
plot("ll", f, cbw, "-o")
// Plot data from figure 12
plot(f0, cbw0,'r')
title("Critical bandwidth vs. frequency (Zwicker et al., 1957)")
legend("Data from table I", "Data from Fig, 12", 2);
xlabel("f  (Hz)")
ylabel("cbw  (Hz)")
gca().grid = [1 1];


scf(2);
clf(2);
// Plot data from figure 12
plot("ll", f0, cbw0, "r")
title("Critical bandwidth vs. frequency (Zwicker et al., 1957)")
legend("Data from Fig, 12", 2);
xlabel("f  (Hz)")
ylabel("cbw  (Hz)")
gca().grid = [1 1];

// Approximation with a least-square spline. It is a 
// spline with given breakpoints (in general as few as
// possible, and in general,  less than the available
// data)
x = linspace(log10(f(1)), log10(f($)), 6)
[y, d] = lsq_splin(log10(f), log10(cbw), x);

fp = f(1) * 10.^linspace(0,log10(f($)/f(1)), 20);
xp = log10(fp);
yp = interp(xp, x, y, d);
cbw_sp = 10.^yp;


scf(3);
clf(3);
plot("ll", fp, cbw_sp,'b')
plot("ll", f0, cbw0, "r")
title("Critical bandwidth vs. frequency (Zwicker et al., 1957)")
legend("LSQ Spline approximation to table I", "Data from Fig, 12", 2);
xlabel("f  (Hz)")
ylabel("cbw  (Hz)")
gca().grid = [1 1];


scf(4);
clf(4);
// Plot standard data from Zwicker
plot('ll', f_zw, cbw_zw, "-o")
// Plot data from the figure 12
plot(f0, cbw0,'r')
title("Critical bandwidth - Comparison between 1961 table and 1957 figure")
xlabel("f  (Hz)")
ylabel("cbw   (Hz)")
legend("1961 standard data","1957 data from figure 12", 2);


// Approximation of the Zwicker 1961 standard data with 
// a least-square spline. It is a spline with given 
// breakpoints (in general as few as possible, and in
// general, less than the available data)
x_zw = linspace(log10(f_zw(1)), log10(f_zw($)), 3);
[y_zw, d_zw] = lsq_splin(log10(f_zw), log10(cbw_zw), x_zw);

fp_zw = f_zw(1) * 10.^linspace(0,log10(f_zw($)/f_zw(1)), 20);
xp_zw = log10(fp_zw);
yp_zw = interp(xp_zw, x_zw, y_zw, d_zw);
cbw_zw_sp = 10.^yp_zw;


scf(5);
clf(5);
plot("ll", f_zw, cbw_zw, "o")
plot("ll", fp_zw, cbw_zw_sp, "r")
title("Approximation of Zwicker 1961 standard data with least-square spline")
gca().data_bounds = [10 50; 13500 4000];
xlabel("f  (Hz)")
ylabel("cbw   (Hz)")
legend("1961 standard data","Least-square spline fitting", 2);
