
// Filtro para ecualizar un auricular audiométrico TDH39.
//
// El filtro está especificado como un circuito RLC pasivo
// en "Psychoacoustics. Facts and models" de Zwicker y 
// Fastl, pág 7. 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. Se cumple Z1*Z2 = R^2, 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.; Fastl, H. (1999) "Psychoacoustics, facts and
//      models". Springer 
//
// --------------------------------------------
// Federico Miyara
// 2020-03-13
// 2024-08-10

clear

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

// Datos del circuito
R = 600;
L1 = [2.65 0.0318];
L2 = [0.00669 0.00795];
C1 = [18.4e-9 22.1e-9];
C2 = [7.3e-6 88.4e-9];
R1 = [10070 2774];
R2 = [36 130];

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 sección 1
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:length(N)
    z(k,:) = roots(N(k))';
    p(k,:) = roots(D(k))';
end


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

// Respuesta en frecuencia compleja de las secciones
for k=1:length(N)
    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 TDH39")
legend("Respuesta completa", "Sección 1", "Sección 2", 4)
xlabel("f   (Hz)")
ylabel("|H|   (dB)") 
gca().grid = [1 1];

