%  This script solves the lecture note example of deriving
%  a model of a room impulse response from the data using Proneys
%  method (IIR).
%
%   Written by Kevin D. Donohue (March 2004).



%  Prony model parameters:
numord = 22;  %  Numerator order
denord = 26;  %  Denomenator order


%  Read in impulse response recording to examine
%  and determine what data segement to present to the algorithm
[y,fs]  =wavread('offroomimp2.wav');
figure(1)
plot(y)
xlabel('Samples')
ylabel('Amplitude')
title('Extract a segment corresponding to single impuse response')
pause
%  Extract impulse response segment
yt = y(58120:60000);
%yt = [zeros(2700-2670,1); y(2700:4000)];
%  Plot segment specta and filter out noise
figure(2)
[p,f] = pwelch(y,hamming(256),128,1024,fs);
plot(f,10*log10(fs*p))
xlabel('Hz')
ylabel('dB')
title('Data Spectrum')
pause
%  From graph oveserved impulse sound from other noise
%  and filter out noise with a bandpass butterworth filter
[b,a] = butter(6, [100, 2500]/(fs/2));
yt = filtfilt(b,a,yt);
figure(3)
plot(yt)
xlabel('Samples')
ylabel('Amplitude')
title('filtered impuse response')
pause
%  Apply prony to impuse segment and plot results

[bp,ap] = prony(yt,numord,denord);
g = filter(bp,ap,[1, zeros(1,length(yt))]);
figure(4);
plot(yt/max(abs(yt)),'r'); hold; plot(g/max(abs(g)),'b');  hold off
title('Impulse response compared (data in red, model in blue)')



%  Examine magnitude and group delay spectra of filter

[h,fa] = freqz(bp,ap,1024,fs);
figure(5); plot(fa,abs(h))
title('Filter Spectra')
[gh,gfa] = grpdelay(bp,ap,1024,fs);

figure(6); plot(gfa,abs(gh))
title('Group Delay')