% BRIR
close all
[y,fs]=audioread('BRIR0.wav');
y=y;%+randn(size(y))*.001;
[b,a]=butter(4,4000*[2^-(1/6) 2^(1/6)]/(fs/2));
Nfft=2^ceil(log2(length(y)));
H=fft(b(:),Nfft)./fft(a(:),Nfft);
y = ifft(fft(y(:,1),Nfft).*H);

% plot(db(y.^2)/2)
hold on
edc=flipud(db(cumsum(flipud(y.^2))))/2;
plot(edc)
% plot(flipud((cumsum(flipud(y.^2))))/2)
% plot((db(cumsum((y.^2))))/2)
ylim([-70 10])

n=3000:30000;
% y = k * n + d
% y = [n(:) ones(N,1)]*[k;d]

X = [n(:), ones(length(n),1)];
kd=X \ edc(n,1);
hold on
plot(n,X*kd,'k--')
text(4000,0,sprintf('RT=%4.1fs',-60/kd(1)/fs))