
// Prueba de la función iso3382.sci
//
// --------------------------------
// Federico Miyara
// 2024-08-14

mode(0)
clear

// Tasa de muestreo
Fs = 44100
// Rango dinámico útil
d = 30
// Fracción de octava
k = 1
// Duración
TT = 2.5
// Vector de tiempos
to = [0:TT*Fs]/Fs;
// Tiempo de reverberación
To = 2
// Constante de tiempo correspondiente al tiempo
// de reverberación
tau = To / log(1000)
// Respuesta al impulso según un modelo simple
// de caída exponencial independiente de la
// frecuencia más un ruido de fondo
ximp = exp(-to/tau).*rand(to, "normal") + 0.001*rand(to, "normal");
// Muestra de ruido de fondo
SN = 0.001*rand(1,44100, "normal");

// Carga de la función iso3382.sci.
exec iso3382.sci;

// Cálculo del tiempo de reverberación total
// sin filtrar por bandas
[Tc, foc, Xc, Pc, Ec, Ic] = iso3382(ximp, Fs, 0, d);

// Vector de tiempos
t = (0:length(Ec)-1)/Fs;

// Error relativo porcentual
er = (Tc/To - 1)*100

scf(1);
clf(1);
gcf().figure_position = [220, 130];
plot(t, 10*log10(Ec))
title("Integración inversa de la respuesta al impulso completa");
xlabel("t  (s)");
ylabel("E  (dB)");
gca().grid = [1 1];

scf(2);
clf(2); 
gcf().figure_position = [240, 110];
plot(t,              10*log10(Ec), ...  
     t(Ic(1):Ic(2)), Pc(1) + Pc(2)*(t(Ic(1):Ic(2)) - t(Ic(1))), "r", ...
     t(Ic(1)),       Pc(1) + Pc(2)*(t(Ic(1)) - t(Ic(1))), "or", ...
     t(Ic(2)),       Pc(1) + Pc(2)*(t(Ic(2)) - t(Ic(1))), "or")
title("Integración inversa de la respuesta al impulsocompleta y ajuste");
xlabel("t  (s)");
ylabel("E  (dB)");
legend("Integración inversa", "Ajuste lineal");
gca().grid = [1 1];


// Cálculo del tiempo de reverberación por bandas.
// Los argumentos de salida también incluyen
// el vector de frecuencias centrales de las 
// bandas, la matriz de las respuestas al impulso
// por bandas, la matriz de parámetros, la matriz
// de las caídas energéticas obr¿tenidas por 
// integración inversa y la matriz de los índices
// donde cada caída energética alcanza los -5 dB y
// los -35 dB respecto a los respectivos máximos
tic
[T, fo, X, P, E, I] = iso3382(ximp, Fs, k, d);
toc

// Errores relativos porcentuales
er = (T/To - 1)*100;

[fo, T, er]

// Graficas de la respuesta al impulso filtrada 
// por bandas de octava
scf(3);
clf(3); 
gcf().figure_position = [260, 90];
for h=1:length(fo)
    subplot(3,3,h)
    plot(t(32768:$), X(h,32768:$))
    legend("fo = " + string(fo(h)));
end

// Graficas de las caídas energéticas de la 
// respuesta al impulso filtrada por bandas 
// de octava
scf(4);
clf(4); 
gcf().figure_position = [280, 70];
for h=1:length(fo)
    t1 = t(I(h,1));
    t2 = t(I(h,2));
    p1 = P(h,1);
    p2 = P(h,2);
    subplot(3,3,h)
    plot(t(32768:$),       10*log10(E(h,32768:$)), ...
    t(I(h,1):I(h,2)), p1 + p2*(t(I(h,1):I(h,2)) - t1), "r", ...
    t(I(h,1)),        p1 + p2*(t(I(h,1)) - t1), "or", ...
    t(I(h,2)),        p1 + p2*(t(I(h,2)) - t1), "or")
    legend("fo = "+string(fo(h)), 3);
end

