% GEN_HT generates a whitening impulse response for a given colored noise
% input file. This file is supposed to be recorded by a non-band limited
% microphone from a white noise source.
%
% Applying the generated impulse response to the data will then yield a
% straight frequency response.
%
% version: 0.1
% 20051018, Remco Stoutjesdijk
% version: 0.2
% 20051019: converted to stereo, verbose via waitbar, interactive file
% selections
%
% TODO:
% - anti-alias
% - indicate error to user
% - compile to c-code / executable

function gen_ht();

close all; clear;

% User interaction for file locations and number of coefficients
[FileName,PathName] = uigetfile('*.wav','Select the input wave file');
if isequal(FileName,0)
   errordlg('No input file selected, aborting')
   return;
end
infile          = strcat(PathName, FileName);

[FileName,PathName] = uiputfile('impulsresponse.wav','Save output file name');
if isequal(FileName,0)
   errordlg('No output file selected, aborting')
   return;
end
outfile         = strcat(PathName, FileName);

answer = inputdlg('Please enter the number of coefficients desired', 'Determine resolution', 1, cellstr('1000'));
N_filter = str2num(char(answer));

% input data
N               = wavread(infile, 'size');
[y, fs, nbits]  = wavread(infile);

if N(2) > 1
   mesg = (['Input file contains ', num2str(N(2)), ' channels.']);
else
   mesg = (['Input file contains ', num2str(N(2)), ' channel.']);
end
h = waitbar(1, mesg); pause(1);

% autocorrelation
waitbar(0, h, 'Calculating autocorrelation');
for channel = 1 : N(2)
   rho(:, channel) = xcorr(y(:, channel));
   waitbar(channel/N(2), h);
end
% figure(1); plot(rho); hold on

% create whitening filter
waitbar(0, h, 'Calculating coefficients');
% N_filter        = 1e3;
for channel = 1 : N(2)
   [coeffs(:, channel), error(channel)] = levinson(rho(N(1):length(rho), channel), N_filter);
   waitbar(channel/N(2), h);
end

nbits           = 32;

mesg = (['Writing wave file: ', num2str(fs), ' Hz, ', num2str(nbits), ' bits, ', num2str(N_filter), ' coeffients...']);
waitbar(0, h, mesg);
% normalize to prevent clipping
coeffs          = coeffs ./ max(max(coeffs));
wavwrite(coeffs, fs, nbits, outfile);
pause(1); close(h);

msgbox(['Whitening impulse response of ' num2str(N_filter) ' coefficients written to ' outfile '.'], 'Done');