% John Nolting
% 424
% Spring Reverb T60

clear all;

[out fs] = wavread('test_sweep_2nd_half.wav');
[out2 fs2] = wavread('spring_2_rch_7_2nd_half.wav');


c = out(1:end);
c_len = length(c);

r = out2(1:end);
% option for uncut waveforms from cubase:
% r = out2(1048575:(1048575+1048574));

r_len = length(r);
ir = real(ifft(fft(r)./fft(c)));

plot(ir)
figure;
minPlotDB = -60;
Bn = [];
An = [];
h_db = [];
%f_list = [125 250 500 1000 2000 4000];
%f_list = [100 500 1000];
f_list = [2000 4000 6000];
for i = 1:length(f_list)
    f = f_list(i)/(fs/2);
    Wn = [0.5*f 2*f];
    [B,A] = BUTTER(2,Wn);
    Bn(:,i) = B;
    An(:,i) = A;
    [h w] = freqz(B,A,1024);
    h_db(:,i) = db(abs(real(h)));
    Y = filter(B, A, ir);
    % Y_len = length(Y);
    % Y_sq = Y*Y;
    
    Y = Y(1:101687);
    nSamp = length(Y);
    

     % calculate the EDC
     signalEDC_DB = zeros(1,nSamp);
     signalEDC_DB(nSamp) = Y(nSamp)^2;
     for n = nSamp-1:-1:1
         signalEDC_DB(n) = signalEDC_DB(n+1) + Y(n)^2;
     end
     signalEDC_DB = 20*log10(signalEDC_DB);


     % normalize the EDC to 0 dB
     maxSignalEDC_DB = max(signalEDC_DB);
     signalEDC_DB_N = signalEDC_DB - maxSignalEDC_DB;
    
     % plot the energy decay curve
     timeAxis = (0:length(Y)-1)/fs; % time axis in seconds
     subplot(3,1,i)
     plot(timeAxis,signalEDC_DB_N);grid;
     % v = axis; axis([v(1) v(2) minPlotDB v(4)]); % to reduce vertical axis
     axis([-inf inf -60 0]);
     title_s = sprintf('Band-passed Impulse from Sweep , Plot to Find T60 Freq = %d Hz',f_list(i));
     title(title_s)
     % title('Normalized Energy Decay Curve (EDC)');
     xlabel('time (s)');ylabel('magnitude (dB)');
     % zoom on;
     orient tall;
end





