% deconvolve.m
%
% Given a system impulse response h[n] and test response y[n], this
% algorithm deconvolves the test response with the impulse response
% to determine the input x[n] as X(jw) = Y(jw)/H(jw).
%
% Eliminate some noise, but is still very rough
%
% NOTE: Assumes all signals are the same length. Use sigcorrect.m
% to equalize their lengths if they are different.
%
% use: deconvolve(impresp,test,testresp)
% variables: impresp: impulse response of system
%            test: test input to system
%            testresp: test response of system
%
% William Howison
% 12/11/2005

function deconvolve(impresp,test,testresp)

N = length(real(fft(impresp))); % determines the length of the signals

% creates a filter to roughly eliminate noise based on previous
% knowledge and tests on our test inputs
firstlast = ones(5e4,1);
middle = zeros(N-2*length(firstlast),1);

filter2 = [firstlast; middle; firstlast];
sinclength = (2*pi)/(N-1);
x = [-pi:sinclength:pi];
filter = fftshift(sinc(x).^2)';
% selects output location to write the deconvolved sound file to
disp('Please select output file location')
[out_file, out_path] = uiputfile('*.wav', 'Select Output File');

HR = fft(impresp); % FFT of the impulse response of the system

T = fft(test); % FFT of original input
TR = fft(testresp); % FFT of the response of the system to the test signal

decon = (TR ./ HR); % naive deconvolution by division of FFTs

%size(decon)
%size(filter)

%decon = decon.*filter; % rough filtering using above created filter
decon = decon.*filter2;

% iFFT the deconvolved signal to get an approximation of the input
decontime = real(ifft(decon));

% record the result (the approximated input) as a wave file
wavwrite(decontime,44100,[out_path out_file]);

% plot inputs and results for comparison versus time
timelength = N/44100; % time length of signals
increment = 1/44100; % increment to measure time axis by
time = increment:increment:timelength; % create time axis

%subplot(4,1,1)
plot(time, test) % plots the original input recorded
title('Original input')

%subplot(4,1,2)
figure(2)
plot(time, decontime) % plots the deconvolved result
title('Deconvolved result (approx. input)')

%subplot(4,1,3)
figure(3)
semilogy(real(abs(T))) % plot real part of original input's spectrum
title('Real spectrum magnitude of original input')

%subplot(4,1,4)
figure(4)
semilogy(real(abs(decon))) % plot real part of result's spectrum
title('Real spectrum magnitude of deconvolved result')