clear
close all

addpath(genpath('Data'));
addpath(genpath('Functions'));
addpath(genpath('Objects'));

load garazs_1.mat

%%
fs = 1 / mean(diff(t));
Nt = 5;
h = h(1:Nt*fs,:);
%
N = 2;
rmic = 0.1;
micArray = spherical_design('File','firthenmike.dat',rmic);
Psi_mic = h(:,1);

T = 0.0025;
Krms = T*fs;    
alpha = exp(-1/Krms);
RMS = sqrt(filtfilt(1/Krms, [1, -alpha], abs(Psi_mic(:,1)).^2));
[~, ix_onset] = max(abs(RMS));
ir_in = Psi_mic(ix_onset:end,1);
t = (0:length(ir_in)-1)'/fs;
freq = (0:length(ir_in)-1)'/length(ir_in)*fs;
EDC_full = flipud(cumsum(flipud(ir_in.^2),1));
EDC_full = 10*log10(EDC_full/EDC_full(1));

ixs = find(EDC_full<=-5 & EDC_full>=-25);
p_full = polyfit(t(ixs),EDC_full(ixs),1);

%
Nc = 25;
fc = 40*2.^((0:Nc)/3);
fl = fc/2.^(1/6);
fu = fc*2.^(1/6);
ir_filt = zeros(length(ir_in),Nc + 1);
for n = 1 : Nc + 1
    [B,A] = butter(2,[fl(n),fu(n)]/(fs/2));
    ir_filt(:,n) = filtfilt(B,A,ir_in);
end
EDC = flipud(cumsum(flipud(abs((ir_filt)).^2),1));
EDC = 10*log10(EDC./EDC(1,:));
figure;
t1 = (0:size(ir_filt,1)-1)/fs;
subplot(2,1,1)
plot(t1, abs((ir_filt(:,1:5:15))))
xlim([0,0.5])

subplot(2,1,2)
plot(t1, EDC(:,1:5:15))
xlim([0,0.5])
%%
 for n = 1 : Nc + 1

    ixs2 = find(EDC(:,n)<=0 & EDC(:,n)>=-10);
    [ p,S] = polyfit(t(ixs2),EDC(ixs2,n),1);
    m(n) = p(1);
    b(n) = p(2);
    RT20oct(n) = -60/m(n);
end

%%
Nstep = 512;
Nhann = 512;
win0 = hann(Nhann);
ix_start = floor(linspace(1,floor(size(ir_in,1)-Nhann-1),Nstep)');
waterfall = zeros(size(ir_in,1),Nstep);
win_mx = zeros(size(ir_in,1),Nstep);
for n = 1 : Nstep
    n

    win = [zeros(ix_start(n),1); 
            win0(1:Nhann/2); 
            ones(size(ir_in,1)-Nhann-ix_start(n),1) ;
            win0(Nhann/2+1:Nhann)];
    win_mx(:,n) = win;
    waterfall(:,n) = 20*log10(abs(fft(ir_in.*win)));
end
%%


%%
fig =  figure('units','normalized','outerposition',[0 0 1 1]);
pos1 = [ 0 0.55 0.5 0.4 ];
pos2 = [ 0.5 0.7 0.5 0.3 ];
pos3 = [ 0.5 0.37 0.5 0.3 ];
pos4 = [ 0.5 0.05 0.5 0.3 ];
pos5 = [0,0.1, 0.5,0.4];


a1 = axes('Units','normalized','outerPosition',pos1);
surf(fc, t(1:1000:end), EDC(1:1000:end,:));view([-250,45]);
xlabel('f_c [Hz]')
ylabel('t [s]');
h=gca;
set(h,'xscale','log')
title('EDC(f_c,t)')

a2 = axes('Units','normalized','outerPosition',pos2);
plot(t,EDC_full,'LineWidth',1.2);
hold on
plot(t,p_full(1)*t+p_full(2),'LineWidth',1.2)
grid on
xlabel('t [s]')
ylabel('EDC [dB]')
yl = ylim;
plot(t([ixs(1),ixs(1)]),yl,'--k')
plot(t([ixs(end),ixs(end)]),yl,'--k')
xlim([t(1),t(end)])
ylim([-100,0])
RT20 = -60/p_full(1);
title(sprintf('Estimated T_{60} = %.3f sec',RT20))

a3 = axes('Units','normalized','outerPosition',pos3);
semilogx(fc,RT20oct,'LineWidth',1.2);
grid on
xlim([fc(1),fc(end)])
xlabel('f [Hz]')
ylabel('T_{60} [s]')
title('Octave-band rev. time')

a4 = axes('Units','normalized','outerPosition',pos4);
H = abs(fft(ir_in));
semilogx(freq,20*log10(H./max(H)),'LineWidth',1.2);
grid on
xlim([fc(1),fc(end)])
xlabel('f [Hz]')
ylabel('T_{60} [s]')
ylim([-40,0])
title('RIR spectrum')

q = 100;
f0 = freq(1:q:end/2,:);
t0 = ix_start/fs;
w2 = waterfall(1:q:end/2,:);
a5 = axes('Units','normalized','outerPosition',pos5);
surf(f0,t0,w2.');
view([-250,45]);
zlim([-90,0])
h=gca;
set(h,'xscale','log')
title('Waterfall diagram')
xlabel('f [Hz]')
ylabel('T [s]')
shading flat

