function [d,sr,SNR] = otanalyze(varargin)
% [d,sr,SNR] = otanalyze(varargin)
% Utility to analyze noisy audio signals signals
%
% Usage;
%   otanalyze -mix mixfile -clean cleanfile
%         -noiseout noisefile -filterout filterfile -targetout targetfile 
%            .. returns separated noise & filter & prints SNR
%
%   See http://labrosa.ee.columbia.edu/projects/otanalyze
%   for full documentation
%
% 2011-08-21 Dan Ellis dpwe@ee.columbia.edu
% $Header: $

VERSION = 0.6;
DATE = 20120120;

% Parse out the optional arguments
[CLEANIN, ...
 TARGETOUT, ...
 MIXIN, ...
 NOISEOUT, ...
 FILTEROUT, ...
 START, END, ...
 SAMPLERATE, ...
 CLEANSKIP, ...
 DISP, ...
 XTRA] = ...
    process_options(varargin, ...
                    '-clean', '', ...
		    '-targetout', '', ...
                    '-mix', '', ...
                    '-noiseout', '', ...
                    '-filterout', '', ...
		    '-start', 0, '-end', 0, ...
                    '-samplerate', 0, ...
                    '-cleanskip', 0, ...
		    '-disp', 0);

HELP = 0;
if length(XTRA) > 0
  HELP = length(strmatch('-help',XTRA,'exact')) > 0;
  if ~HELP
    disp(['Unrecognized options:',sprintf(' %s',XTRA{1:end})]);
    HELP = 1;
  end
end

  
if HELP
  disp(['otanalyze v',num2str(VERSION),' of ',num2str(DATE)]);
  disp('usage: otanalyze ...');
  disp('   -clean <filename>  The name of the clean (reference) sound file');
  disp('   -targetout <filename>  Where to save extracted target speech');
  disp('   -mix <filename>    Input mixture of noise and channel-filtered target');
  disp('   -noiseout <filename> Where to write the separated noise residual');
  disp('   -filterout <filename> Where to save the estimated channel response');
  disp('   -start <time_secs> Start processing at this time in the files');
  disp('   -end <time_secs>   Finish processing at this point in the files');
  disp('   -cleanskip <time_secs>  Drop time from start of CLEAN (or MIX if -ve)');
  disp('   -samplerate <sr_hz> Resample to (low) sampling rate before xcorr');
  disp('   -disp <bool>       display graphics if set');

  if nargout > 1
    d = []; sr = []; SNR = [];
  end

  return
end

if SAMPLERATE == 0
  sr = []; 
else
  sr = SAMPLERATE;
end

MONO = 1; % mix down before processing

disp(['++++++ otanalyze v',num2str(VERSION),' ++++++']);

%ndisp = 5;
ndisp = 4;
dispn = 0;

% CLEANSKIP +ve means skip sound at start of CLEAN
if CLEANSKIP > 0
  SKIPC = CLEANSKIP;
  SKIPM = 0;
else
  SKIPC = 0;
  SKIPM = -CLEANSKIP;
end

% mix + clean to give noise + filter + SNR
disp(['Reading MIX from ',MIXIN,' ...']);
DUR = max(0,END-START);
[dmix,sr] = audioread(MIXIN, sr, MONO, START+SKIPM, DUR);

%wavwrite(dmix,sr,'dmix.wav');

disp(['Reading CLEAN from ',CLEANIN,' ...']);
[dclean,sr] = audioread(CLEANIN, sr, MONO, START+SKIPC, DUR);

%wavwrite(dclean,sr,'dclean.wav');

disp('Identifying CLEAN in MIX...');

% first-pass xcorr to find mix in clean
[n,xc] = find_skew(dclean, dmix);
if n > 0
  % mix starts midway through clean
  % chop off first portion of clean
  % so that clean starts just 1/2 sec before mix comes in
  % and ends 1/2 sec after end of mix 
  cleanmin = max(0, n-round(sr/2));
  cleanlen = min(length(dclean)-cleanmin, length(dmix)+sr);
  dclean = dclean(cleanmin + [1:cleanlen]);
else
  % mix starts before clean
  % chop off first part of mix
  % so that first part of mix is 1/2 sec into clean
  % and mix lasts no longer than clean
  mixmin = (-n)+round(sr/2);
  cleanmin = -mixmin;
  mixlen = min(length(dmix)-mixmin,length(dclean));
  dmix = dmix(mixmin + [1:mixlen]);
end

[dnoise, dtarg, dfilt, SNR, delay] = ...
    find_in_mix(dmix, dclean, sr);
% add on what we took off
tdelay = delay + (cleanmin/sr) + CLEANSKIP;
disp(['Mix delay= ',sprintf('%.3f',tdelay),' s']);
%disp(['  = ',sprintf('%.3f',-delay),' + ',sprintf('%.3f',cleanmin/sr), ...
%      ' + ',sprintf('%.3f',CLEANSKIP)]);

dispn = maybedisp(dmix,sr,dispn,ndisp,DISP,...
                  ['Mix: delay=',sprintf('%.3f',tdelay), ' s, '...
                   'SNR=',sprintf('%.1f',SNR),' dB - ',MIXIN]);

disp(['Mix SNR= ',sprintf('%.2f',SNR),' dB']);

if nargout > 0
  d = dnoise;
end
  
dispn = maybedisp(dnoise,sr,dispn,ndisp,DISP,['noise - ',NOISEOUT]);
dispn = maybedisp(dtarg,sr,dispn,ndisp,DISP,['target - ',TARGETOUT]);

newgain = 1;
if ~isempty(FILTEROUT)
  if max(abs(dfilt)) > 1
    newgain = 32767 / 32768 / max(abs(dfilt));
    disp(['Warning: reducing gain of filter (and noise, and target) by ', ...
          num2str(newgain),' x to avoid clipping.']);
  end
  audiowrite(newgain*dfilt,sr, FILTEROUT);
  disp(['FILTER saved to ',FILTEROUT]);
end
if ~isempty(NOISEOUT)
  audiowrite(newgain*dnoise, sr, NOISEOUT);
  disp(['NOISE saved to ',NOISEOUT]);
end
if ~isempty(TARGETOUT)
  audiowrite(newgain*dtarg,sr, TARGETOUT);
  disp(['TARGET saved to ',TARGETOUT]);
end

if DISP
  if dispn <= ndisp
    dispn = dispn+1;
    subplot(ndisp,2,2*dispn-1);
    plot([1:length(dfilt)]/sr,dfilt);
    title(['filter - ',FILTEROUT]);
    [H,W] = freqz(dfilt);
    HdB = 20*log10(abs(H));
    subplot(ndisp,2,2*dispn);
    plot(W/pi*sr/2,HdB);
    axis([0, sr/2, max(HdB) + [-60 0]]);
    grid;
    title(['filter freq response']);
  end
end

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  dispnout = maybedisp(d,sr,dispnin,ndisp,DISP,tit)
% Maybe add a spectrogram and move down
%fmax = 4000;
dbrange = 60;

% 32 ms FFT window len
fftwinsec = 0.032;
fftlen = 2^(round(log(sr*fftwinsec)/log(2)));

dispnout = 0;

if DISP
  if dispnin < ndisp
    dispnout = dispnin+1;
    subplot(ndisp,1,dispnout);
    specgram(d,fftlen,sr);
%    ax = axis();
%    ax(4) = fmax;
%    axis(ax);
    ca = caxis();
    ca(1) = ca(2)-dbrange;
    caxis(ca);
    colorbar;
    title(tit,'interpreter','none');
  end
end