function [decay_curve, EDT, T60] = IR_Analyze(x,fs)
 
% IR_Analyze
%
% Impulse Response Analyzer -
%
%   Reads and analyzes the audio file of an impulse response. Produces
%   a representative backwards-integrated Schroeder Decay Curve, EDT & T60
%   measurements
%
%   Input Parameters:
%
%       x  = mono audio signal/vector
%       fs = sampling frequency
%
%   Ouput Parameters:
%
%       decay_curve = column vector representing normalized decay values
%                     of input signal, in decibels
%       EDT         = Early time decay: (time of signal decay from 0dB to
%                     -10 dB) * 6
%       T60         = 2*[(time of -60dB) - (time of -5dB)]
%    
 
 
%======================================================
% Determine number of channels within audio file; reduce to mono if the
% number of channels exceed 1
[len,channels] = size(x);
if channels > 1
    z = zeros(len,1);
    for i = 1:channels-1
        z = z + x(:,1);
    end  
    x = z;
end
 
%======================================================
% Generate backwards-integrated Schroeder Decay Curve
[max_x n0] = max(x); 
curve = cumsum(x(end:-1:n0).^2);
curve = curve(end:-1:1);
 
% Convert curve values to decibel values
curve_db = 20*log10(abs(curve));
% Normalize decay curve to zero
norm_curve_db = curve_db - max(curve_db);
decay_curve = norm_curve_db;
 
 
%======================================================
% Plot Schroeder Decay Curve
figure; 
timeAxis = (0:length(curve)-1)/fs;
plot(timeAxis,decay_curve);grid
 
% Label X and Y axis
ylabel('dB');
xlabel('time(s)');
 
% Determine boundaries for  X and Y axis
ylim([-70 0]);                    % Y axis ends at -70 dB
 
xMax = find(decay_curve < -60);   % X axis ends where decay curve
xMax = min(xMax)/fs;              % meets -60dB
xlim([0 xMax]);
 
%======================================================
% Determine EDT
edt = find(decay_curve > -10);
EDT = (max(edt) / fs) * 6;
% Determine T60
t60 = find(decay_curve > -60);
t60 = (max(t60))/fs;
t5db = find(decay_curve > -5);
t5db = (max(t5db))/fs;
T60 = 2*(t60 - t5db);