
// Parametrización de la curva de altura psicoacústica
// (ratio pitch) de Zwicker-Fastl
// 
// Se parte de la curva proveniente de la figura 5.5 de 
// "Psychoacoustics. Facts and models" de Zwicker-Fastl
// (1999, p. 112). Representa en eje doble logarítmico
// zm (altura psicoacústica) vs. f (frecuencia) entre
// los límites 10 Hz - 20000 Hz y 20 mel - 2469 mel. La 
// curva debe pasar por (fo, zmo) = (125 Hz, 125 mel).
// La gráfica fue copiada al portapapeles con la 
// herramienta snapshot de Acrobat y pegada en Word. Un
// problema es que la imagen original escaneada tiene 
// una leve rotación de 0.29° en sentido antihorario, 
// por lo que hubo que rotarla para compensar. Para 
// ello se pegó como imagen nueva en GIMP y se rotó con
// Image / Transform / Arbitrary rotation, se guardó y  
// abrió en Word, se calcó con la herramienta de dibujo 
// de Word, se separó de la imagen original, se pasó
// a PDF, se abrió en Acrobat, con Snapshot se copió
// la curva al portapapeles y se abrió en GIMP, donde
// se recortó para ajustar los límites para que la 
// esquina inferior izquierda y la superior derecha 
// coincidieran con los extremos de la curva. Luego se
// pasó a modo Greyscale y se exportó como bitmap. 
//
// Finalmente se compara con el índice de banda crítica 
//
//
// --------------------------------
// Federico Miyara
// 2024-10-15
// 2024-10-28

clear

// Ruta al archivo que contiene la gráfica
path = "R:\Mis documentos\Acústica\Ruido\" + ..
       "PSICOACÚSTICA COMPUTACIONAL\_GRÁFICAS\"
path = "R:\Mis documentos\Acústica\Ruido\" + ..
       "PSICOACÚSTICA COMPUTACIONAL\PSICOACÚSTICA\"
arch = "Ratio_pitch_x_log20_log20000_y_log20_log2469.bmp"
bmpfile = path + arch;

// Límites de la grafica ubicando el origen de
// coordenadas en la esquina inferior izquierda
xmin = log10(20/20)
xmax = log10(20000/20)
ymin = log10(20/20)
ymax = log10(2469/20)

// Cargar función que convierte el mapa de bits
// en tabla
exec bmp2table.sci;

// Obtener la tabla de valores x e y 
tic
[y, x] = bmp2table(bmpfile, xmin, xmax, ymin, ymax);
toc

// Corrección de y(1) ya que la curva real vale
// 0 en la abscisa inicial
y(1) = 0;

// Recuperación de valores de frecuencia y 
// altura psicoacústica sin corregir rotación
f = 20*10.^x;
zm = 20*10.^y;

// Ajustar diferencia debida a resolución
// de la gráfica
[mini, i] = min(abs(f - 125))
zm1 = zm*f(i)/zm(i);

// Cálculo del índice de banda crítica
exec zbark.sci;
z = zbark(f);

// Gráficas
scf(1);
clf(1);
plot("ll", f, zm, f, zm1, 125, 125, "or")
title("Curva de altura psicoacústica en función de la frecuencia")
xlabel("f  (Hz)");
ylabel("zm  (mel)");
legend("zm sin corregir", "zm corregida", "Referencia 125 Hz",2);
gca.data_bounds = [20 20; 20000 3000];
gca.grid = [1 1];
gca().tight_limits = ["on","on","off"];

scf(2);
clf(2);
plot("ln", f, zm, f, z*100, 125, 125, "or")
title("Altura psicoacústica e índice de banda crítica vs. frecuencia")
xlabel("f  (Hz)");
ylabel("zm  (mel)")
legend("zm (mel)", "100*z", "Referencia 125 Hz",2);
h1 = gca();
gca().data_bounds = [20 20; 20000 3000];
gca().tight_limits = ["on","on","off"];
gca().grid = [1 1];
h2 = newaxes();
h2.filled = "off";
h2.axes_visible(2) = "on";
h2.data_bounds = h1.data_bounds
h2.y_location = "right"
ylabel("z  (bark)")
h2.y_ticks.labels = string((0:15)*2)'








