
// Prueba de la función late reverb
//
// En este script se da un ejemplo de aplicación de 
// la función, la aplicación por convolución a una
// señal vocal, la medición mediante el método de
// Schroeder de los tiempos de reverberación junto
// con el cálculo del error relativo para cada 
// frecuencia y las gráficas de  la respuesta al 
// impulso y las caídas energéticas en dB con sus
// respectivos ajustes lineales. 
//
// ----------------------------------------------
// Federico Miyara
// 2015-08-15

clear
mode(0)

// Datos
To = 3
L = [4 6 2.7]
f =     [  63  125  250  500 1000 2000 4000 8000 16000]
alpha = [0.05 0.10 0.12 0.18 0.25 0.30 0.35 0.45  0.50]
Fs = 44100
N = 4096

// Volumen y area de la sala a partir de sus 
// dimensiones
V = L(1)*L(2)*L(3);
S = 2*(L(1)*L(2) + L(1)*L(3) + L(2)*L(3));

// Tiempos de reverberación según Sabine
Tsab = 0.161*V/S./alpha'

// Cálculo de la respuesta al impulso
exec late_reverb.sci;
y = late_reverb(To, L, f, alpha, Fs, N);
t = [0:length(y)-1]/Fs;

// Reproducción de la respuesta al impulso
playsnd(y/10, Fs, 1)

// Carga de un sonido vocal
[v, Fsv] = wavread("G:\Mis documentos\Ondas\hemodinamia.wav");

// Reproducción del sonido
playsnd(v, Fsv, 1)

// Convolución con la respuesta al impúlso
z = convol(v, y);
z = 0.3*z/max(abs(z));

// Reproducción de la señal húmeda
playsnd(z, Fsv, 1)

// Reproduccion de la mezcla
playsnd(z(1:length(v)) + 0.7*v, Fsv)

// Grafica de la respuesta al impulso (sin
// reflexiones tempranas) 
scf(1);
clf(1);
gcf().figure_position = [220 110];
plot(t,y)
xlabel("t  (s)")
ylabel("Impulse response") 

// Cálculo del tiempo de reverberación 
// correspondiente a la respuesta al impulso 
// por el método de Shroeder de la integración
// inversa
exec iso3382.sci;
[T, fo, X, P, E, I] = iso3382(y, Fs, 1, 30);
tE = (0:size(E,2)-1)/Fs;

// Cálculo del error relativo porccntual y tabla
// de valores del tiempo de reverberación de
// Sabine y medido mediante el método de Schroeder
er_perc = (T - Tsab)./Tsab * 100 
disp([Tsab, T, er_perc])

// Respuesta al impulso para las diferentes bandas
// con sus respectivos ajustes lineales
scf(2);
clf(2); 
gcf().figure_position = [240, 90];
for h=1:length(fo)
    t1 = tE(I(h,1));
    t2 = tE(I(h,2));
    p1 = P(h,1);
    p2 = P(h,2);
    subplot(3,3,h)
    plot(tE(32768:$),       10*log10(E(h,32768:$)), ...
    tE(I(h,1):I(h,2)), p1 + p2*(tE(I(h,1):I(h,2)) - t1), "r", ...
    tE(I(h,1)),        p1 + p2*(tE(I(h,1)) - t1), "or", ...
    tE(I(h,2)),        p1 + p2*(tE(I(h,2)) - t1), "or")
    legend("fo = "+string(fo(h)));
end

