function s060719_timedomain_chnuss( par )

% function s060719_timedomain_chnuss( par )
%
% Mandatory parameters:
%   par.input_filename
%   par.input_format
%   par.do_chn
%   par.do_uss
%   par.output_filename
%
% Action:
% - Read a waveform "par.input_filename" in the format
%   "par.input_format" ('wav', 'sphere', 'raw-le', 'raw-be').
% - Filter it, applying CHN (iff "par.do_chn == 1") followed by USS
%   (iff "par.do_uss == 1").
% - Write out a .wav file "par.output_filename".
%
% There are many other optional parameters, including
% "par.frameshift_sec" and "par.framelength_sec" (default 0.010 and
% 0.025). See the code for their description:
%    dbtype s060719_timedomain_chnuss
% 
% See IDIAP RR 06-09.
%
% By Guillaume Lathoud - lathoud@idiap.ch - 2006.
  
  if nargin < 1
    error( [ mfilename ' needs 1 input arguments.' ] );
  end
  
  required_fields = { 'input_filename', 'input_format', 'do_chn', 'do_uss', 'output_filename' };
  check_param( required_fields, fieldnames( par ) );

  % Check for input file
  if ~exist( par.input_filename, 'file' )
    error( [ mfilename ': input file does not exists ("' par.input_filename '").' ] );
  end

  % Check for output file
  if exist( par.output_filename, 'file' )
    disp( [ mfilename ': output file already exists ("' par.output_filename '"). Skipped.' ] );
    return;
  end
  
  % Check the input format

  par.input_format = lower( par.input_format );
  
  if ~ismember( par.input_format, { 'wav', 'sphere', 'raw-le', 'raw-be' } )
    error( [ mfilename ': unknown input format.' ] );
  end
  
  % Get sampling frequency
  
  if ismember( par.input_format, { 'raw-le', 'raw-be' } )
    % From the user
    check_param( { 'fs' }, fieldnames( par ) );
  end

  if strcmp( par.input_format, 'wav' )
    % From the file itself
    [ a_size, par.fs ] = wavread( par.input_filename, 'size' );
  end

  if strcmp( par.input_format, 'sphere' )
    % From the file itself
    sphere_header = read_sphere_header( par.input_filename );
    par.fs = sphere_header.sample_rate;
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 2 ) Fill in default parameters
  
  %%%%%%%%%%
  % Main parameters
  
  par_default = [];
  par_default.frameshift_sec  = 0.010;
  par_default.framelength_sec = 0.025;
  
  % Size of the block on which we apply CHN / USS
  par_default.block_size_sec = 1.0;
  
  % Debugging: show windows and do "keyboard" pauses
  par_default.debug = 0;
  
  % Verbosity
  par_default.verbose = 0;
  
  %%%%%%%%%%
  % Secondary parameters
  par_default = fill_second( par_default );
  
  
  par = fill_default( par, par_default );


  % Default verbosity for RSE fitting
  if isempty( par.rse_em_par.verbose )
    par.rse_em_par.verbose = max( 0, par.verbose - 1 );
  end
  
  % Deduce values in samples
  par.framelength = round( par.framelength_sec * par.fs );
  par.frameshift  = round( par.frameshift_sec  * par.fs );
  
  % From "par.framelength_sec", deduce the block size in number of frames
  
  par.block_size_frames = max( 1, round( par.block_size_sec / par.frameshift_sec ) );
  
  % Parameters for "reduced" data representation
  
  if isempty( par.future_nframes_max )
    par.future_nframes_max = par.block_size_frames;
  end
  
  if isempty( par.past_nframes_max )
    par.past_nframes_max   = par.future_nframes_max;
  end


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 3 ) Open input and output files
  
  % Input: open and get number of samples
  % For "raw-le", "raw-be" and "wav", 16-bit samples are assumed
  
  fid_in = [];
  
  switch par.input_format
    
   case 'raw-le'
    fid_in = fopen( par.input_filename, 'r', 'ieee-le' );
    fseek( fid_in, 0, 'eof' );
    par.nsamples = floor( ftell( fid_in ) / 2 );
    fseek( fid_in, 0, 'bof' );
    
   case 'raw-be'
    fid_in = fopen( par.input_filename, 'r', 'ieee-be' );
    fseek( fid_in, 0, 'eof' );
    par.nsamples = floor( ftell( fid_in ) / 2 );
    fseek( fid_in, 0, 'bof' );
    
   case 'wav'
    a_size = wavread( par.input_filename, 'size' );
    par.nsamples = a_size( 1 );
    
   case 'sphere'

    sphere_header = read_sphere_header( par.input_filename );
    par.nsamples  = sphere_header.sample_count;

    switch sphere_header.sample_byte_format

     case '01'
      fid_in = fopen( par.input_filename, 'r', 'ieee-le' );

     case '10'
      fid_in = fopen( par.input_filename, 'r', 'ieee-be' );

     otherwise
      error( [ mfilename ': unknown sphere_header.sample_byte_format "' sphere_header.sample_byte_format '"' ] );

    end

   otherwise
    error( [ mfilename ': input format "' par.input_format '" not supported yet.' ] );
  end
  
  % From "par.nsamples", deduce the number of frames
  
  par.nframes = 0;

  if par.nsamples >= par.framelength
    
    par.nframes = 1 + floor( ( par.nsamples - par.framelength ) / par.frameshift );
    
  end
  
  % From "par.fs", deduce the band on which we fit the RSE model
  
  par.rse_band = max( 1, min( par.framelength, round( par.rse_band_hz / par.fs * 2 * par.framelength ) ) );

  % Open output file

  fid = fopen( par.output_filename, 'wt' );
  fclose( fid );

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 4 ) Read the input signal
  
  switch par.input_format
    
   case { 'raw-le', 'raw-be' }
    
    x_in = fread( fid_in, par.nsamples, 'int16' );
    
   case 'wav'
    
    x_in = wavread( par.input_filename, par.nsamples );
    
   case 'sphere'
    
    fseek( fid_in, sphere_header.header_size, 'bof' );
    x_in = fread( fid_in, a_block.nsamples, 'int16' );
    
   otherwise
    
    error( [ mfilename ': input format "' par.input_format '" not supported yet.' ] );
    
  end
  
  fclose( fid_in );

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 5 ) Filter the signal
  
  x_out = chnuss_signal_main( x_in, par );
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 6 ) Write out the filtered signal

  wavwrite( x_out / ( eps + 1.01 * max( abs( x_out ) ) ), par.fs, 16, par.output_filename );

  if par.verbose
    disp( [ mfilename ' wrote out "' par.output_filename '".' ] );
  end
  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define default secondary parameters

function par_default = fill_second( par_default )

  % Pre-emphasis
  par_default.preemp = 0.97;  
  
  % Zero-mean activated or not
  par_default.do_zero_mean = 1;
  
  %%%%%%%%%%
  % Parameters for data representation
  par_default.past_nframes_max   = [];
  par_default.future_nframes_max = [];
  
  %%%%%%%%%%
  % Parameters for channel estimation

  par_default.chn_alpha = 0.20;  % Use the 20% lowest magnitude values

  %%%%%%%%%%
  % Parameters for fitting the Rayleigh + Shifted Erlang model

  par_default.rse_band_hz = [ 300  3400 ];  % Quite useful on Aurora
  
  par_default.rse_em_par.erlang_t_factor = 1;
  par_default.rse_em_par.min_iter        = 5;
  par_default.rse_em_par.max_iter        = 200;
  par_default.rse_em_par.sigma_floor     = 1e-10;
  par_default.rse_em_par.lambda_floor    = 1e-10;
  
  % 0: moment method (faster and better)
  % 1: ML method (exact)
  par_default.rse_em_par.do_fminsearch = 0;
  
  opt = optimset( @fminsearch );
  opt = optimset( opt, 'Display', 'none', 'MaxFunEvals', 50, 'MaxIter', 25, 'TolX', 1e-4, 'TolFun', 1e-4 );
  par_default.rse_em_par.fminsearch_options = opt;    

  % Activate default verbosity
  par_default.rse_em_par.verbose = [];
