
// Análisis de un modelo circuital de la cóclea
//
// El modelo consiste en una línea de transmisión con 
// parámetros variables con la distanoa a la ventana
// oval
//
// Referencia
// Deutsch, Sid; Micheli-Tzanakou, Evangelina (1987)
//    "Neuroelectric Systems". New york University
//    Press. New York 
// Deutsch, Sid; Deutsch, Alice (1993) "Understanding 
//    the nervous system. An engineeering perspective".
//    IEEE Press.New York
// Miyara, Federico (2024) Análisis de un modelo circuital
//    de la cóclea 
//
// ----------------------------------------
// Federico Miyara
// 2024-09-26
// 2024-09-30

clear
close(winsid());

// Carga de la función cochlea_model()
exec cochlea_model.sci;

// Carga de la función cochlea_travelling_wave()
exec cochlea_travelling_wave.sci;

// ----------------------------------------
// 1 ENVOLVENTES
// ----------------------------------------

// Vector de distancias a la ventana oval
x = 0:0.0002:0.035;
caso = 2
select caso
case 1
   // Vector de frecuencias utilizadas en Deutsch
   // et al. (1994)
   f = [170 1160 4700 15000]';
case 2
   // Otro vector de frecuencias
   f = [200 1000 5000 15000]';
end

// Obtener las respuestas espaciales para
// cada frecuencia 
[H, f] = cochlea_model(x, f);

// Obtener los máximos de cada envolvente
// y los índices donde éstos se dan
[Hmax, imax] = max(abs(H),"c")

// Posiciones tonotópicas de los máximos
xmax = x(imax)

// Normalización de las envolventes respecto
// a sus respectivos máximos
Hnor = H ./ repmat(Hmax,1,length(x));

// Gráfica de las envolventes
scf(1);
clf(1);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
gcf().figure_size = [800,550];
// Título contextualizado
tit1 = "Envolventes para ";
for i=length(f):-1:2
    tit1 = tit1 + string(f(i)) + " Hz, ";
end    
tit1 = tit1 + string(f(1)) + " Hz";
title(tit1)
plot(x, abs(Hnor))
legend(string(f'));
xlabel("x  (m)")
ylabel("y_env")
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, 0; 0.040, 1.1];


// ----------------------------------------
// 2 ONDA VIAJERA
// ----------------------------------------

// Vector de distancias a la ventana oval
x = 0:0.00001:0.035;

// Frecuencia
f = 1000;

// Fase
alpha = 0;

// Vector de tiempos
T = 1/f;
t = 0:T/6:5/6*T;

// Obtención de la onda viajera en varios
// instantes sucesivos y la envolvente
[y, env] = cochlea_travelling_wave(x, t, f, alpha);

// Gráfica de la deformación de la membrana
// basilar a intervalos de 1/6 de ciclo
scf(2);
clf(2);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
gcf().figure_size = [800,550];
plot(x, env, "k--", x, -env, "k--", x, y)
tit2 = "Patrones de oscilación de la membrana basilar para ";
tit2 = tit2 + string(f) + " Hz en instantes sucesivos";
title(tit2)
xlabel("x  (m)")
ylabel("yBM")
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, -1.1; 0.035, 1.1];

// Serie de gráficas conteniendo una curva por
// cada gráfica
scf(3);
clf(3);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15 ];
gcf().figure_size = [850,550]
for i=1:6
   subplot(2,3,i)
   plot(x, env, "k--", x, -env, "k--", x, y(i,:))
   if i==1
      lbl = "t = 0";
   elseif i==2
      lbl = "t = T/3";
   else
      lbl = "t = " + string(i-1) + "T/3";
   end
   xstring(0.027, 0.7, lbl, 1);
   xlabel("x (m)")
   gca().data_bounds = [0, -1.1; 0.035, 1.1];
   gca().tight_limits = ["on","on","on"];
end


// ----------------------------------------
// 3 VELOCIDAD DE LA ONDA VIAJERA
// ----------------------------------------

// Vector de distancias a la ventana oval
deltax = 0.0001;
x = 0:deltax:0.035;

// Frecuencias
f = [100, 200, 1000, 2000]; 

// Fases
alfa = ones(f)*%pi/4;

// Dos tiempos que difieren en un pequeño incremento
deltat = 0.0001;
t = [0, deltat]

// Obtención de la onda viajera en dos instantes
// próximos para dos frecuencias diferentes
[y, env] = cochlea_travelling_wave(x, t, f, alfa); 

// Figura a exportar para ilustrar el concepto
// de velocidad de fase  
scf(4);
clf(4);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15 ];
gcf().figure_size = [800,550];
plot(x, env(3,:), "k--", x, -env(3,:), "k--", x, y(:,:,3))
tit2 = "Patrones de oscilación de la membrana basilar para ";
tit2 = tit2 + string(f(3)) + " Hz en instantes sucesivos";
title(tit2)
xlabel("x  (m)")
ylabel("yBM")
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, -1.1; 0.035, 1.1];

// Se toman ahora gran cantidad de puntos 
// en x para tener buena precisión en la 
// derivada 
deltax = 0.00001;
x = 0:deltax:0.035;

// Función de transferencia de la cóclea
// para las frecuencias requeridas
H = cochlea_model(x, f);

// Envolventes
env = abs(H) ./ (max(abs(H),"c") * ones(1,size(H,2)));

// Fase de la función de transferencia en 
// radianes
phi = phasemag(H, "c")*%pi/180;

// Cálculo de la velocidad de fase de la onda
// viajera. Aplicaremos la fórmula 
// c = -2*%pi*f / (dfase/dx)
c = -2*%pi*repmat(f(:), 1, length(x)-1) ./ (diff(phi,1,"c")/deltax);
xc = x(1:length(x)-1);

scf(5);
clf(5);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
gcf().figure_size = [650,550];
subplot(2,1,1)
plot(xc, c)
title("Velocidad de fase de la onda viajera en la cóclea")
legend(string(f) + " Hz");
//xlabel("x  (m)")
ylabel("c  (m/s)")
subplot(2,1,2)
plot(x, env)
xlabel("x  (m)")
ylabel("y_env")
gca().tight_limits = ["on","on","on"];
gca().data_bounds = [0, 0; 0.035, 1.1];
 
// ----------------------------------------
// 4 MÁXIMOS
// ----------------------------------------

// Longitud de la membrana basilar
L = 0.035
// Resistencia característica de la cóclea
Ro = 0.02
// Coeficiente de los exponentes de w
a1 = -100
a2 = -20
// Coeficientes de las partes real e imaginaria
// de la impedancia por unidad de longitud
A1 = 0.39
A2 = 0.1325
// Coeficientes de las partes real e imaginaria
// del exponente de la corriente del modelo  
B1 = A1/Ro*2/a1^2
B2 = A2/Ro*2/a2^2

// Posición
xx = 0.002:0.0005:0.035;

// Transformación no lineal de la posición
w = sqrt(L - xx);  

// Coeficientes de la ecuación cuadrática de
// la frecuencia al cuadrado en función de la
// transformación no lineal de la posición 
P = A1^2*2*B1*a1^2*exp(2*(a1-a2)*w).*exp(a1*w).*w;
Q = 2*a1*A1^2*exp(2*(a1-a2)*w) + A2^2*2*B1*a1^2*exp(a1*w).*w;
R = 2*a2*A2^2;

// Resolución de la ecuación cuadrática
f = sqrt((-Q + sqrt(Q.^2 - 4*P*R))/2./P);

// Curva tonotópica obtenida por mí a partir 
// de la curva de Deutsch et al. (1993)
ff  =  10.^(1.3012919 + 0.3863057*(35.034435 - xx*1000).^0.5652075);

// Error
er = abs(f - ff);
// Eliminación del error infinito para x = L
er_rel = er(1:$-1)./f(1:$-1);
// Error relativo máximo
[er_rel_max, i] = max(er_rel);
f_er_rel_max = f(i);
disp(f_er_rel_max, er_rel_max)

// Gráfica en función de la frecuencia
scf(6);
clf(6);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
gcf().figure_size = [650,550];
subplot(2,1,1)
plot("nl", xx, f, xx, ff)
title("Frecuencia para la cual el máximo se da en x") 
legend("f_modelo", "f_tonotópica")
xlabel("xmax  (m)");
ylabel("f  (Hz)"); 
gca().grid = [1 1];
gca().data_bounds = [0, 20; 0.035, 20000];
gca().tight_limits = ["on","on","on"];
subplot(2,1,2)
plot (xx(1:$-1), 100*er_rel)
title("Error entre el modelo y la curva tonotópica real") 
gca().grid = [1 1];
xlabel("x  (m)");
ylabel("error relativo  (%)");
gca().data_bounds = [0, 0; 0.035, 100];
gca().tight_limits = ["on","on","on"];


// ----------------------------------------
// 5 RESONANCIA EN UNA POSICIÓN ESPECÍFICA
// ----------------------------------------

// Aproximación de la posición del máximo para 
// fo = 1000 
fo = 1000;
[delta, i] = min(abs(f - fo));
f(i-1:i+1)
xx(i-1:i+1)
xo = interp1(f(i-1:i+1), xx(i-1:i+1), 1000, "spline")

fff = 20*10.^(0:0.01:3);  
H = cochlea_model(xo, fff);
HdB = 20*log10(abs(H));
[Hmax, i] = max(HdB)
disp(fff(i), Hmax)

scf(7);
clf(7);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
plot("ln", fff, HdB, fff(i), Hmax, "or")
title("Resonancia en la posición de máxima oscilación para " + .. 
      string(fo) + " Hz")
xlabel("f  (Hz)");
ylabel("|H(f)|  (dB)");
gca().grid = [1 1];
gca().data_bounds = [20, max(HdB)-80; 20000, max(HdB)+10];
gca().tight_limits = ["on","on","on"];


// ----------------------------------------
// 6 MÚLTIPLES TONOS
// ----------------------------------------

// Aquí combinaremos múltiples tonos. Supondremos
// que las resonancias respecto a la posición se
// encuentran reforzadas a través de un exponente e
e = 8;

// Vector de distancias a la ventana oval
x = 0:0.00001:0.035;

// Frecuencias
f = [400 1000 2500];

// Fases
alpha = zeros(f);

// Obtención de la respuesta espacial compleja
HH = cochlea_model(x, f);

// Envolventes normalizadas
env = abs(HH)./ (max(abs(HH),"c") * ones(1,size(HH,2)));

// Envolvente global original
env_global = sqrt(sum(env.^2,"r"));

// Aplicación de una transformación no lineal
// que emula una acentuación del pico de 
// resonancia
env1 = env.^8;

// Envolvente global transformada
env_global1 = sqrt(sum(env1.^2,"r"));

scf(8);
clf(8);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15 ];
subplot(2,1,1)
plot(x, abs(env1), x, env_global1, "k--")
title("Envolventes con picos acentuados y su envolvente global")
xlabel("x  (m)")
ylabel("y_env")
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, 0; 0.035, 1.6];
subplot(2,1,2)
plot(x, abs(env), x, env_global, "k--")
title("Envolventes con los picos originales y su envolvente global")
xlabel("x  (m)")
ylabel("y_env")
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, 0; 0.035, 1.6];

// Fase en radianes
phi = phasemag(HH, "c")*%pi/180;

// Vector de tiempos
t = 0:0.00005:0.005;

// Obtener las ondas viajeras para todos los
// instantes y todas las frecuencias
y = cochlea_travelling_wave(x, t, f, alpha, e); 

// Sumar las ondas de las diversas frecuencias
yy = sum(y,3);

// Graficar las ondas viajeras para todos los 
// instantes de tiempo en que se las calculó
mini = min(yy)-0.1;
maxi = max(yy+0.1);
scf(9);
clf(9);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
subplot(2,1,1)
plot(x, yy)
title("Ondas viajeras para diversos instantes")
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, mini; 0.035, maxi];
xlabel("x  (m)")
ylabel("yBM")
subplot(2,1,2)
plot(x, env_global1, "r", x, -env_global1, "r") 
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, mini; 0.035, maxi];
xlabel("x  (m)")
ylabel("env_global(x)")

scf(10);
clf(10);
fig = gcf().figure_id;
gcf().figure_position = [180 + fig*20, 160 - fig*15];
plot(x, yy(1,:))
plot(x, env_global1, "r--", x, -env_global1, "r--") 
gca().tight_limits = ["on","on","off"];
gca().data_bounds = [0, -3; 0.035, 3];
