// Function iso3382.sci test
//
// --------------------------------
// Federico Miyara
// 2024-08-14

mode(0)
clear

// Sample rate
Fs = 44100
// Useful dynamic range
d = 30
// Octava fraction
k = 1
// Duration
TT = 2.5
// Time vector
to = [0:TT*Fs]/Fs;
// Reverberation time
To = 2
// Time constant corresponding to reverberation
// time
tau = To / log(1000)
// Impulse response from a simple exponential decay
// model independent of frequency with a background
// noise
ximp = exp(-to/tau).*rand(to, "normal") + 0.001*rand(to, "normal");
// Background noise sample
SN = 0.001*rand(1,44100, "normal");

// Load function iso3382.sci
exec iso3382.sci;

// Calculate total reverberation time without any
// band filtering
[Tc, foc, Xc, Pc, Ec, Ic] = iso3382(ximp, Fs, 0, d);

// Time vector
t = (0:length(Ec)-1)/Fs;

// Percent relative error 
er = (Tc/To - 1)*100

// Unfiltered impulse response energy decay 
scf(1);
clf(1);
gcf().figure_position = [220, 130];
plot(t, 10*log10(Ec))
title("Reverse integration of the full-band impulse response");
xlabel("t  (s)");
ylabel("E  (dB)");
gca().grid = [1 1];

// Impulse response with linear fitting
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("Reverse integration of the impulse response and its fitting");
xlabel("t  (s)");
ylabel("E  (dB)");
legend("Reverse integrafion", "Linear fitting");
gca().grid = [1 1];


// Band reverberation time
// Output arguments also include centre frequency
// vector, band impulse response matrix, band 
// parameter matrix, band energy decay matrix and
// matrix including, for each band, the indices of
// the instants when the energy decay falls 5 dB 
// and 35 dB below its maximum  
tic
[T, fo, X, P, E, I] = iso3382(ximp, Fs, k, d);
toc

// Percent relative errores 
er = (T/To - 1)*100;

[fo, T, er]

// Octave-band filtered impulse response plots
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

// Energy decay of the octave-band-filtered  
// impulse response 
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

