% John Nolting
% 424
% Spring reverb
% Uses butter function to create a set of butterworth filters.  
% The filters are octave-wide, bandpass, and fourth order 
% Calculates snr and returns to plot title

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);
% r = out2(1048575:(1048575+1048574));

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

% plot(ir)

Bn = [];
An = [];
h_db = [];
f_list = [100 500 1000 2000 4000 6000];
%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);
    inits = (sqrt(sum(Y(600:800,1).^2)/201));
    % noise = db(0);
    %pulse_height = db(sqrt(sum(Y(1:15).^2)/(15)));
    noise = db(sqrt(sum(Y(488800:498800,1).^2)/10001));
    inits    
    pulse_height = db(max(Y));
    snr = pulse_height - noise;
    snr
    %subplot(3,1,i)
    %figure;
    plot(Y)
    title_s = sprintf('Impulse Response Passed Thru Butterworth with Freq = %d Hz, SNR = %d dB',f_list(i), round(snr));
    title(title_s)
    orient tall;
end

% plot butter bands
figure;
semilogx(w, h_db(:,1), w, h_db(:,2), w, h_db(:,3), w, h_db(:,4), w, h_db(:,5), w, h_db(:,6));
axis([0 3 -60 0]);
title('Transfer function Magnitudes of Butter Band Passes')
xlabel('Samples (Log Scale)')
ylabel('Amplitude (dB)')