// Cálculo de las curvas de igual nivel de sonoridad
// por los métodos de la IRAM 4066 : 1970 (ISO 226 
// anterior a 1970), ISO 226 : 1987 e ISO 226 : 1997(?) 

clear, clc, close(winsid())

// ---------------------------------------------
// ISO/R 226:1961 - IRAM 4066:1970
// Corresponde a las curvas de Robinson-Dadson
// Proporciona una ecuación cuadrática para 
// calcular el nivel de sonoridad en función 
// del nivel de presión sonora
// ---------------------------------------------

// Datos correspondientes a IRAM 4066:1970
// Ver al final
f70 = [     20     25     30     35    40    45    50    55    60 ..
            70     80     90    100   120   140   160   180   200 ..
           250    300    350    400   500   600   700   800   900 ..
          1000   1200   1400   1600  1800  2000  2500  3000  3500 ..
          4000   4500   5000   6000  7000  8000  9000 10000 11000 ..
         12000  13000  14000  15000];
af70 = [-217.2 -167.2 -135.7 -113.2 -96.3 -83.2 -73.0 -65.1 -58.7 ..
         -49.1  -42.5  -37.4  -33.5 -27.0 -22.7 -19.4 -16.9 -14.7 ..
         -10.8   -8.1   -6.1   -4.7  -3.0  -1.8  -1.0  -0.5  -0.2 ..
           0.0    0.5    1.1    1.7   2.5   3.3   5.3   6.9   7.9 ..
           7.9    7.1    5.3   -0.5  -7.5 -13.3 -16.5 -16.8 -14.8 ..
         -12.7  -13.9  -22.7  -43.0 ] 
bf70 = [ 3.669  3.145  2.801  2.538 2.325 2.154 2.021 1.918 1.836 ..
         1.723  1.652  1.603  1.570 1.512 1.473 1.444 1.422 1.404 ..
         1.362  1.325  1.290  1.259 1.205 1.155 1.109 1.064 1.028 ..
         1.000  0.972  0.953  0.937 0.927 0.924 0.928 0.937 0.946 ..
         0.954  0.963  0.973  1.011 1.075 1.159 1.242 1.314 1.377 ..
         1.450  1.566  1.777  2.146];
cf70 = [  -927   -782   -677   -589  -515  -454  -407  -370  -341 ..
          -305   -285   -274   -269  -260  -253  -248  -244  -242 ..
          -231   -218   -202   -185  -151  -117   -84   -50   -22 ..
             0     23     47     71    90   100   118   120   113 ..
            98     80     59     14   -35   -89  -145  -203  -269 ..
          -350   -454   -591   -772]/1e5;
Luf70 = [ 74.3   65.0   58.1   52.7  48.4  44.8  41.7  39.1  36.8 ..
          32.9   29.8   27.2   25.1  21.4  18.9  16.8  15.2  13.8 ..
          11.2    9.4    8.1    7.2   6.0   5.2   4.7   4.4   4.3 ..
           4.2    3.9    3.4    2.6   1.8   1.0  -1.2  -2.9  -3.9 ..
          -3.9   -3.0   -1.1    4.6  10.9  15.3  17.0  16.4  14.2 ..
          12.0   12.0   16.0   24.1];

// Frecuencias equirrepartidas logarítmicamente
ff70 = 20*10^(0:0.05:log10(15000/20));

// Cálculo de 12 curvas de igual nivel de sonoridad
// Nivel de presión sonora en función de la frecuencia
// para cada nivel de sonoridad entre 10 fon y 120 fon
for k=1:12
   for i=1:length(f70)
      // El modelo matemático da el nivel de sonoridad
      // en función del nivel de presión sonora como un
      // polinomio de segundo grado, por lo que para 
      // obtener el nivel de presión sonora a partir de
      // el nivel de sonoridad es necesario despejar 
      // el nivel de presión sonora correspondiente a 
      // un nivel de sonoridad mediante la resolvente 
      // de una ecuación de segundo grado. 
      A = cf70(i);
      B = bf70(i);
      C = af70(i) - k*10;
      Lpf70(k,i) = (-B + sqrt(B^2 - 4*A*C))/2/A;
      // Corregir el valor en 1000 Hz, donde la
      // resolvente no es aplicable pues A=0 
      if f70(i)==1000
         Lpf70(k,i) = k*10;
      end
   end 
end

// Cálculo de 12 curvas de igual nivel de presión 
// sonora. Nivel de sonoridad en función de la 
// frecuencia para cada nivel de presión sonora 
// entre 10 dB y 120 dB
for k=1:12
   for i=1:length(f70)
      Lnf70(k,i) = af70(i) + bf70(i)*10*k + cf70(i)*(10*k)^2;
   end
end

scf(1);
clf(1);
gcf().figure_position = [200,0]
gcf().figure_size = [837,735]
subplot(2,2,1)
// Lp, parámetro Ln
plot("ln", f70, Lpf70)
plot("ln", f70, Luf70, "--r")
title("Curvas isofónicas según IRAM 4066:1970")
xlabel("f  (Hz)")
ylabel("Lp  (dB)")
gca().data_bounds = [20 -10; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
subplot(2,2,3)
// Ln, parpametro Lp
plot("ln", f70, Lnf70)
title("Curvas de igual nivel de presión sonora según IRAM 4066:1970")
xlabel("f  (Hz)")
ylabel("Ln  (fon)")
gca().data_bounds = [20 0; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
// Parámetros af, bf, cf
subplot(3,2,2)
plot("ln", f70, af70)
title("Parámetro a(f)")
xlabel("f  (Hz)")
gca().grid = [1 1];
gca().data_bounds = [20 -220; 20000 20];
gca().tight_limits = ["on","on","on"];
subplot(3,2,4)
plot("ln", f70, bf70)
title("Parámetro b(f)")
xlabel("f  (Hz)")
gca().grid = [1 1];
gca().data_bounds = [20 0; 20000 4];
gca().tight_limits = ["on","on","on"];
subplot(3,2,6)
plot("ln", f70, cf70)
title("Parámetro c(f)")
xlabel("f  (Hz)")
gca().grid = [1 1];
gca().data_bounds = [20 -0.01; 20000 0.002];
gca().tight_limits = ["on","on","on"];




// ---------------------------------------------
// ISO 226:1987 - IRAM 4066:1997
// Corresponde a las curvas de Robinson-Dadson
// Proporciona una ecuación para calcular el 
// nivel de presión sonora en función del nivel 
// de sonoridad
// ---------------------------------------------

// Datos corresondientes a ISO 226:1987 
// Frecuencias centrales normalizadas de tercio de octava
f87  = [   20    25  31.5    40    50    63    80   100   125    160 ..
          200   250   315   400   500   630   800  1000  1250   1600 ..
         2000  2500  3150  4000  5000  6300  8000 10000 12500];
af87 = [2.347 2.190 2.050 1.879 1.724 1.597 1.512 1.466 1.426  1.394 ..
        1.372 1.344 1.304 1.256 1.203 1.135 1.062 1.000 0.967  0.943 .. 
        0.932 0.933 0.937 0.952 0.974 1.027 1.135 1.266 1.501];
bf87 = [ 5.61  5.27  4.81  4.04  3.38  2.86  2.59  2.57  2.56   2.55 ..
         2.54  2.48  2.29  2.01  1.62  1.11  0.52  0.00 -0.39, -0.67 ..
        -0.92 -1.05 -1.04 -0.88 -0.55  0.00  0.89  2.11  4.88] * 0.001;
Tf87 = [ 74.3  65.0  56.3  48.4  41.7  35.5  29.8  25.1  20.7   16.8 ..
         13.8  11.2   8.9   7.2   6.0   5.0   4.4   4.2   3.7    2.6 ..
          1.0  -1.2  -3.6  -3.9  -1.1   6.6  15.3  16.4  11.6]; 

// Con este juego de parámetros se puede calcular 
// el nivel de sonoridad Ln a partir del nivel de 
// presión sonora Lpf a la frecuencia f así:
//
//               af(Lpf - Tf)
// Ln = 4.2 + -----------------
//             1 + bf(Lpf - Tf)  
//
// o bien puede calcularse el nivel de presión sonora 
// correspondiente a un nivel de sonoridad mediante
//
//                Ln - 4.2
// Lpf = Tf + -------------------
//            af - bf(Ln - 4.2)

// Cálculo de 51 frecuencias logarítmicamente
// equiespaciadas entre 20 Hz y 12,5 kHz
ff87 = 20 * 10.^(log10(12500/20)*[0:0.02:1]);

// Interpolación de los parámetros af, bf, Tf
aff87 = interp1(f87, af87, ff87, "spline");
bff87 = interp1(f87, bf87, ff87, "spline");
Tff87 = interp1(f87, Tf87, ff87, "spline");

// Se calculan 11 curvas
for k=1:11
   // Un cálculo para cada frecuencia
   Ln = k*10;  
   for i=1:length(ff87)
      Lpf87(k, i) = Tff87(i) + ..
                    (Ln - 4.2) / (aff87(i) - bff87(i)*(Ln - 4.2));
   end
end

// Gráficas
scf(2);
clf(2);
subplot(1,2,1)
// Gráficas de Lp para 11 valores de Ln desde
// 10 dB hasta 110 dB 
plot("ln", ff87, Lpf87)
// Gráfica del umbral
plot("ln", ff87, Tff87, "--r")
title("Curvas de igual nivel de sonoridad ISO 226:1987")
xlabel("f  [Hz]")
ylabel("Lp  [dB]")   
gca().grid = [1 1];
gca().data_bounds = [20 -10; 20000 120];
gca().tight_limits = ["on","on","on"];
// Gráficas de los parámetros af, bf, Tf interpolados
subplot(3,2,2)
plot("ln", ff87, aff87)
title("Parámetro af - ISO 226:1987")
xlabel("f   [Hz]")
ylabel("af")   
gca().grid = [1 1];
gca().data_bounds = [20 0; 20000 2.5];
gca().tight_limits = ["on","on","on"];
subplot(3,2,4)
plot("ln", ff87, bff87)
title("Parámetro bf - ISO 226:1987")
xlabel("f   [Hz]")
ylabel("bf")   
gca().grid = [1 1];
gca().data_bounds = [20 -0.002; 20000 0.006];
gca().tight_limits = ["on","on","on"];
subplot(3,2,6)
plot("ln", ff87, Tff87)
title("Parámetro Tf - ISO 226:1987")
xlabel("f   [Hz]")
ylabel("Tf   (dB)")   
gca().grid = [1 1];
gca().data_bounds = [20 -10; 20000 80];
gca().tight_limits = ["on","on","on"];


// ------------------------------------------------
// ISO 226:2003 - IRAM 4066:2008
// Corresponde a las curvas de Suzuki-Takeshima (?)
// Proporciona una ecuación para calcular el nivel 
// de presión sonora en función del nivel de 
// sonoridad
// ------------------------------------------------

// Parámetros ISO 226 : 1997 (?)  -  IRAM 4066 : 2004
f03  = [      20    25  31.5    40    50    63    80   100   125 ..
             160   200   250   315   400   500   630   800  1000 ..
            1250  1600  2000  2500  3150  4000  5000  6300  8000 ..
           10000 12500];

alfaf03 = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 ..
           0.330 0.315 0.301 0.288 0.276 0.267 0.259 0.253 0.250 ..
           0.246 0.244 0.243 0.243 0.243 0.242 0.242 0.245 0.254 ..
           0.271 0.301];
LUf03 = [  -31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3  -8.1  -6.2 ..
            -4.5  -3.1  -2.0  -1.1  -0.4   0.0   0.3   0.5   0.0 ..
            -2.7  -4.1  -1.0   1.7   2.5   1.2  -2.1  -7.1 -11.2 ..
           -10.7  -3.1]; 
Tf03 = [    78.5  68.7  59.5  51.1  44.0  37.5  31.5  26.5  22.1 ..
            17.9  14.4  11.4   8.6   6.2   4.4   3.0   2.2   2.4 ..
             3.5   1.7  -1.3  -4.2  -6.0  -5.4  -1.5   6.0  12.6 ..
            13.9  12.3];

// Con este juego de parámetros se puede calcular
// el nivel de sonoridad Ln a partir del nivel de 
// presión sonora Lf a la frecuencia f así:
//
// Ln = 40*log(Bf) + 94
// 
// donde
//
// Bf = 0.4 * 10.^((Lp + LU -90)/10).^alfaf 
//      - 0.4 * 10.^((TTf + LU -90)/10).^alfaf + 0.005135
//
// o bien se puede calcular el nivel de presión sonora
// correspondiente a un nivel de sonoridad mediante
//
// Lp = (10/alfaf)* log(Af) - LU + 94
//
// donde 
//
// Af = 4.47e-3*(10.^(Ln/40)-1.15) + (0.4*10.^((TTf+LU-90)/10)).^alfaf


// Cálculo de 51 frecuencias logarítmicamente
// equiespaciadas entre 20 Hz y 12,5 kHz
ff03 = 20 * 10.^(log10(12500/20)*[0:0.02:1]);

// Interpolación de los parámetros alfaf, LUf, Tf
alfaff03 = interp1(f03, alfaf03, ff03, "spline");
LUff03 = interp1(f03, LUf03, ff03, "spline");
Tff03 = interp1(f03, Tf03, ff03, "spline");

// Se calculan 10 curvas para niveles de 
// sonoridad entre 10 dB y 100 dB
for k=1:10
   Ln = k*10;
   // Un cálculo para cada frecuencia
   for i=1:length(ff03)
      Af = 4.47e-3 * (10^(0.025*Ln) - 1.15) + ..
           (0.4 * 10^((Tff03(i) + LUff03(i) - 90)/10))^alfaff03(i);
      Lpf03(k,i) = 10/alfaff03(i)*log10(Af) - LUff03(i) + 94;
   end
end


scf(3);
clf(3);
subplot(1,2,1)
// Gráficas para 10 niveles de sonoridad
// entre 10 dB y 100 dB 
plot("ln", ff03, Lpf03)
// Gráfica del umbral
plot("ln", ff03, Tff03, "r--")
title("Curvas de igual nivel de sonoridad ISO 226:2003")
xlabel("f  [Hz]")
ylabel("Lp  [dB]")   
gca().grid = [1 1];
gca().data_bounds = [20 -10; 20000 130];
gca().tight_limits = ["on","on","on"];
// Gráficas de los parámetros
subplot(3,2,2)
plot("ln", ff03, alfaff03)
title("Parámetro alfaf - ISO 226:2003")
xlabel("f   [Hz]")
ylabel("alfaf")   
gca().grid = [1 1];
subplot(3,2,4)
plot("ln", ff03, LUff03)
title("Parámetro LUf - ISO 226:2003")
xlabel("f   [Hz]")
ylabel("alfaf")   
gca().grid = [1 1];
subplot(3,2,6)
plot("ln", ff03 ,Tff03)
title("Parámetro Tf - ISO 226:2003")
xlabel("f   [Hz]")
ylabel("alfaf")   
gca().grid = [1 1];


// ------------------------------------------------
// Aplicación a un ejemplo
// ------------------------------------------------

// Frecuencias centrales de las bandas de tercio de octava
fn  = [   20    25  31.5    40    50    63    80   100   125    160 ..
          200   250   315   400   500   630   800  1000  1250   1600 ..
         2000  2500  3150  4000  5000  6300  8000 10000 12500];

// Índice de banda crítica para las frecuencias centrales
// de las bandas de tercio de octava (según Zwicker-Fastl)
zn = 13 * atan(fn/1316) + 3.5 * atan((fn/7500).^2);

// Cálculo de 100 frecuencias logarítmicamente equiespaciadas
// entre 20 Hz y 12,5 kHz
ff = 20 * 10.^(log10(12500/20)*[0:0.01:1]);

// Índice de banda crítica correspondiente(según Zwicker-Fastl)
zz = 13 * atan(ff/1316) + 3.5 * atan((ff/7500).^2);

// Bandas críticas. Cada columna contiene la frecuencia
// superior y el ancho de banda de cada banda crítica 
bc = [ 100  200   300   400  510  630  770  920 1080 1270 ...
      1480 1720  2000  2320 2700 3150 3700 4400 5300 6400 ...
      7700 9500 12000 15500; ...
       100  100   100   100  110  120  140  150  160  190 ..
       210  240   280   320  380  450  550  700  900 1100 ..
      1300 1800  2500  3500];

// Curva de transferencia para niveles de sonoridad
// Ln <= 40 (no sería necesario, ya que se dispone 
// de la ecuación)
son = [ 8:40; ...
        0.013 0.018 0.023 0.030 0.038 0.046 0.056 0.068 0.080 0.090 ..
        0.105 0.120 0.140 0.157 0.180 0.205 0.235 0.260 0.295 0.330 ..
        0.370 0.410 0.450 0.500 0.540 0.585 0.640 0.690 0.750 0.810 ..
        0.875 0.940 1.000];


// Ejemplo: Onda cuadrada
// Luego puede reemplazarse por cualquier otro espectro,
// donde Lp contendrá los niveles de presión sonora 
// correspondientes a las frecuencias f 
Lptotal = 60;
f1 = 200;
f = [1:2:101] * f1;
Lp = [];
Lp(1) = 60 - 20*log10(%pi*sqrt(2)/4);    
for i=3:2:101
   Lp = [Lp, Lp(1) - 20*log10(i)];
end

scf(4);
clf(4);
// plot as vertical bars
plot2d3(f, Lp)
title('Armónicos de la onda cuadrada')
xlabel('f  [Hz]')
ylabel('Lp  [dB]')   


// Cálculo del Lp en cada banda crítica (LpBC)

n0 = 1;
p2 = zeros(1,24);
for h=1:24
   while f(n0) <= bc(1,h),
      p2(h) = p2(h) + 10^(Lp(n0)/10);   
      n0 = n0+1;
      if n0 > length(f)
         break
      end  
   end
   if n0 > length(f)
      break
   end  
   LpBC(h) = 10*log10(p2(h));
end

disp(LpBC)

// Cálculo de la frecuencia central de las bandas críticas
foBC = bc(1, 1:24) - bc(2, 1:24)/2;

// Cálculo del nivel de sonoridad
// 
//               af(Lf - Tf)
// Ln = 4.2 + -----------------
//             1 + bf(Lf - Tf)  

// Cálculo de la sonoridad por un método simplificado
// no tiene en cuenta la reducción de sonoridad por
// enmascaramiento
Stot = 0;
for h=1:24
   //disp(h)
   if LpBC(h)<0
      Ln(h) = -%inf;
   else
      Ln(h) = interp1(f87, af87, foBC(h), "spline") * ..
              (LpBC(h) - interp1(f87, Tf87, foBC(h), "spline")) / ..
              (1 + interp1(f87, bf87, foBC(h), "spline") * ..
              (LpBC(h) -  interp1(f87, Tf87, foBC(h), "spline"))) + 4.2;
   end
   if Ln(h)<40
      if Ln(h)<2
         S(h) = 0;
      else
         S(h) = 6.901e-5*(Ln(h)-2)^2.634;
      end
   else
      S(h) = 2^((Ln(h) - 40)/10);
   end
   Stot = Stot + S(h);
end

disp([foBC', Ln, S])
disp(Stot)

scf(5);
clf(5);
plot (son(1,:),son(2,:))
title('Curva de transferencia para NS <= 40 fon')
xlabel('NS  [fon]')
ylabel('S  [son]')


// Valores de cf70 aparentemente erróneos en el trabajo original
// de Robinson y Dadson en el que se basa la ISO/R 266:1961 en 
// la que se basa la IRAM 4066:1970. Esos valores fueron 
// cooregidos en la ISO dando curvas más suaves en altos niveles
// de sonoridad y bajas frecuencias 

//cf70 = [  -927   -782   -677   -589  -514  -454  -408  -371  -343 ..
//          -305   -285   -275   -268  -260  -255  -249  -245  -241 ..
//          -230   -218   -204   -188  -151  -116   -84   -48   -22 ..
//             0     23     47     71    90   100   118   120   113 ..
//            98     80     59     14   -35   -89  -145  -203  -269 ..
//          -350   -454   -591   -772]/1e5;
