% deconv_impulse.m

% Takes a recording of impulse response signal and
% "deconvolves" the signal by dividing the corresponding
% FFTs.  The FFT of impulse response is multiplied by a 
% complex constant to correct for time lag in the impulse
% response recording.

% NOTE: assumes all files same sample length to ensure that FFTs
% are the same length too!

clear all;

Fs = 44100;             % sampling rate of audio data
N = 705600;             % length of each audio sample
cutoff = 88100;         % cutoff sample for lowpass filter
lowpass = [ones(cutoff,1); zeros(N-(2*cutoff),1); ones(cutoff,1)];


disp('loading audio data...');

[h_file, h_path] = uigetfile('*.wav', 'Select Impulse Response');
[test_file, test_path] = uigetfile('*.wav', 'Select Audio Data');
[out_file, out_path] = uiputfile('*.wav', 'Select Output File');

h = wavread([h_path h_file],N);           % import impulse response data
test = wavread([test_path test_file],N);  % import test audio                                        

% perform phase correction for impulse response time delay
disp('normalizing and phase correcting impulse response...');

h = h ./ (pi*max(abs(h(:,2))));  % normalize impulse response signal
H = fft(h);
H0(:,1) = H(:,1) .* conj(H(:,2));
H0(:,2) = H(:,2) .* conj(H(:,2));

h0 = real(ifft(H0(:,1)));

% deconvolve impulse response from zero padded test data
disp('deconvolving impulse response from test audio...');
result = real(ifft((fft(test(:,1)) ./ H0(:,1)) .* lowpass));

disp('making plots...');

figure(1);

time = 16/N : 16/N : 16;
subplot(4,1,1);
plot(time, h0);
title('Dabney Lounge Impulse Response');
subplot(4,1,2);
plot(time, test(:,2));
title('Input Signal');
subplot(4,1,3);
plot(time, test(:,1));
title('Recorded Output Signal');
subplot(4,1,4);
plot(time, result);
title('Deconvolved Output Signal');

figure(2);
semilogy(abs(fft(result)));
title('|FFT| of output signal');

disp('playing deconvolved output signal...');
soundsc(result, Fs);

disp('writing results to output file...');
wavwrite(result,Fs,[out_path out_file]);


disp('clearing used variables...');
clear result lowpass h h0 test time;
clear N cutoff H H0;

disp('done.');