
// Filtro para ecualizar un auricular audiométrico DT48.
//
// El filtro está especificado como un circuito RLC pasivo
// en "The ear as a cmmunication receiver" de Zwicker y 
// Feldtkeller, pág 51. Se trata de una red de R constante,
// es decir, formada por secciones tales que la impedancia 
// vista desde la entrada es resistiva pura (independiente
// de la frecuencia) e igual a la resistencia de carga, R. 
// La transferencia de cada sección es Z2/(Z2 + R), donde 
// Z2 es la impedancia de la rama conectada a masa, en este
// caso una serie RLC:
// 
// Z2 = R2 + 1/(C2*s) + L2*s
//
// Resulta, para cada sección 
//
//        R2 + 1/(C2*s) + L2*s          1 + R2*C2*s + L2*C2*s
// H = -------------------------- = -----------------------------
//      R2 + R + 1/(C2*s) + L2*s     1 + (R2 + R)*C2*s + L2*C2*s
//
// Referencia
// Zwicker, E.; Felrtkeller, R. (1999) "The ear as a communication 
//      Receiver". ASA-AIP.
// 
// -------------------------------
// Federico Miyara
// 2020-03-15
// 2024-06-19
// 2024-08-10

clear

// Variables simbólicas (polinomiales)
// s: variable de Laplace
s = %s

// Datos del circuito
R = 600;
L1 = [1.9 0.45 0.0081];
L2 = [0.026 0.007 0.0175];
C1 = [72e-9 20e-9 49e-9];
C2 = [5.2e-6 1.25e-6 22.5e-9];
R1 = [4000 4200 820];
R2 = [90 86 440];

H = 1 + R2.*C2*s + L2.*C2*s^2;
H = H ./ (1 + (R2 + R).*C2*s + L2.*C2*s^2);

// V: señal de entrada
V = 1

// Numerador y denominador de las secciones
N = H.num
D = H.den

// Numerador y denominador normalizado 
N = N/coeff(N)(1)
D = D/coeff(D)(1)

// Ceros y polos de las seeciones
for k=1:3
    z(k,:) = roots(N(k))';
    p(k,:) = roots(D(k))';
end


// Vector de frecuencias
f = 10.^[0:.01:4]*20;
w = 2*%pi*f;

// Respuesta en frecuencia compleja de las secciones
for k=1:3
    HH(k,:) = horner(N(k)/D(k), %i*w);
end

// Respuesta en frecuencia compleja total
HHtot = horner(prod(N./D), %i*w);

// Gráfica
scf(675);
clf(675);
plot("ln", f, 20*log10(abs([HHtot', HH'])))
title("Respuesta del ecualizador para campo libre del DT48")
legend("Total", "Sección 1", "Sección 2", "Sección 3",4)
xlabel("f  [Hz]");
ylabel("|H|  [dB]");
gca().grid = [1 1];
