
// Curva tonotópica según Deutsch & Deutsch
//
// Parametrización y ajuste de la curva de la figura 4.6 de 
// la referencia, que proporciona la frecuencia a la que
// se tiene el pico de respuesta para cada posición x en la
// cóclea 
// 
// Referencias
// Deutsch, Sid; Deutsch, Alice (1994) Understanding tne 
//     nervous system. An Engineering perepective. IEEE Press 
// Zwislocki, J.J. (1965) Analysis of some auditory 
//     characteristics, en andbook of Pathematical Psychology,
//     vol 3, (editores R. D. Luce et al. New York. Wiley.
//
// --------------------------------------
// Federico Miyara
// 2024-08-02

clear

// ------------------------------------------------------
// 1) Datos 
// ------------------------------------------------------

// En el primer intento se parametrizan manualmente 
// algunos puntos de la gráfica de la figura 9.6
//

// Escala de la figura 9.6 respecto a medidas en mm
Cx = 81.5/35;
Cy = 69/log10(20000/20);

// Datos parametrizados de la figura 9.6
x = [ 0.00  5.00 10.00 15.00 20.00 25.00 30.00 32.21 34.03 34.50 35.00]
y = [66.50 60.75 54.50 48.00 40.75 32.25 21.75 16.00  9.25  6.25  0.00]

// Conversión de las ordenadas en frecuencia
f = 20*10^(y/Cy);

// Gráfica de la curva
scf(1);
clf(1);
gcf().figure_position = [210,120];
plot("nl", x, f, x, f, "ro") 
gca().grid = [1 1];
gca().data_bounds = [0 20 ; 35 20000]; 
title("")
xlabel("x  (mm)");
ylabel("f  (Hz)")

// Datos parametrizados automáticamente de la 
// gráfica editada para eliminar grilla
// NOTA: La gráfica debe estar muy limpia, sin
// ruido y sin ningún píxel fuera de la curva.
// El grosor de la curva no es relevante siempre
// que no tenga huecos y no sea exagerado, ya
// que el algoritmo promedia los puntos en cada
// abscisa 
exec bmp2table.sci;
Path = "G:\Mis documentos\";
Path = Path + "Acústica\Ruido\PSICOACÚSTICA COMPUTACIONAL\_GRÁFICAS\";
File = "Resonance_frequency_vs_place_x_0-35_mm_y(log)_20-20000_Hz_low.bmp";
bmpfile = Path + File;
// Límites de la gráfica
xmin = 0;
xmax = 35;
ymin = log10(20);
ymax = log10(20000);
// Obtención de la tabla de datos
tic
[y,x] = bmp2table(bmpfile, xmin, xmax, ymin, ymax);
toc
// Conversión del eje vertical a frecuencia 
f = 10^y;
 
// Gráfica
scf(2);
clf(2);
gcf().figure_position = [230,100];
plot("nl",x,f);
title("Tonotopía en la cóclea: frecuencia pico en función de x")
xlabel("x  (mm)");
ylabel("f   (Hz)");
gca().grid = [1 1]; 

// ------------------------------------------------------
// 2) Aproximación mediante spline por mínimos cuadrados 
// ------------------------------------------------------

// Distancia entre abscisas de los nodos propuestos
// de modo de tener N nodos
N = 21
delta = (max(x) - min(x)) / (N - 1);
// Vector de abscisas de los nodos
xnodo = min(x) + (0:N-1)* delta;
// Obtención de los parámetros de la spline  
[fnodolog, d] = lsq_splin(x(1:10:$), log10(f(1:10:$)), xnodo)
// Valores de abscisa donde se calculará la aproximación
// spline
xs = min(x):0.1:max(x);
// Aplicación de la spline para obtener los valores
// de ordenadas aproximadas
fslog = interp(xs, xnodo, fnodolog, d);
// Conversión a frecuencia
fs = 10^fslog;

// Error
finterps = interp1(x,f,xs);
[maxer, j] = max(abs(fs - finterps))
[maxerr, i] = max(abs(fs - finterps)./fs) 
fs(i)

// Gráfica
scf(3);
clf(3);
gcf().figure_position = [250,80];
subplot(1,2,1) 
plot("nl", x, f, xs, fs) 
gca().grid = [1 1];
gca().data_bounds = [0 20 ; 35 20000]; 
title("Frecuencia pico en función de x")
xlabel("x  (mm)");
ylabel("f  (Hz)")
subplot(1,2,2)
plot("nl", (fs - finterps)./fs*100, fs)
ylabel("f  [Hz]")
xlabel("Error relativo   (%)")
gca().grid = [1 1];


// ------------------------------------------------------
// 3) Ajuste mediante funciónes logarítmica y potencial 
// ------------------------------------------------------
//
// Se utilizarán las funciones 
//
// f(x) = a + b*log10(c - x) 
//
// f(x) = a + b*(c - x)^d 

function e = G(p, z)
    x = z(1,:);
    y = z(2,:);
    e = y - (p(1) + p(2)*log10(p(3) - x));
endfunction

function e = GG(p, z)
    x = z(1,:);
    y = z(2,:);
    e = y - (p(1) + p(2)*abs(p(3) - x).^p(4));
endfunction

// Ajuste logarítmico
Z = [x; log10(f)];
c = 36
a = log10(f($))
b = log10(f(1)/f($))/log10(36)
po = [a, b, c]
[p,err] = datafit(G, Z, po)

fcalc = 10.^(p(1) + p(2)*log10(p(3) - x));

scf(4);
clf(4);
gcf().figure_position = [270,60];
plot("nl", x, f, x, fcalc, "r")
title("Curva tonotópica: ajuste logarítmico")
xlabel("x  (mm)")
ylabel("f   (Hz)")
legend("Datos", "Ajuste logarítmico", 3);
gca().grid = [1 1];


// Ajuste potencial
// Se extrapolan un poco los datos para evitar que
// el punto de pendiente -inf caiga dentro del
// rango de x. Para ello se agrega un x > 35  
xf = 35.05;
yf = log10(f($)) + log10(f($)/f($-1))/(x($)-x($-1))*0.05;
ZZ = [Z, [xf; yf]];
d = 0.6
c = 35.5
b = log10(f(1)/f($))/(c^d - (c-35)^d)
a = log10(f(1)) - b*c^d 
ppo = [a b c d]
[pp,err] = datafit(GG, ZZ, ppo)

ffcalc = 10.^(pp(1) + pp(2)*(pp(3) - x).^pp(4));

scf(5);
clf(5);
gcf().figure_position = [290,40];
plot("nl", x, f, x, ffcalc, "r")
title("Curva tonotópica: ajuste logarítmico")
xlabel("x  (mm)")
ylabel("f   (Hz)")
legend("Datos", "Ajuste potencial", 3);
gca().grid = [1 1];

// Errores del ajuste con función potencial
e = ffcalc - f;
[Emax, imax] = max(abs(e))
fimax = f(imax)

er = (ffcalc - f)./f;
[Ermax, irmax] = max(abs(er))
firmax = f(irmax)

scf(6);
clf(6);
gcf().figure_position = [310,20];
gcf().figure_size = [848,588]
plot("ln", f, er*100)
title("Error relativo en la frecuencia")
xlabel("f  (Hz)")
ylabel("er   (%)")
gca().grid = [1 1];

