function chnuss_x = chnuss_signal_main( x, par )

% function chnuss_x = chnuss_signal_main( x, par )
%
% Time-domain filtering of a noisy speech waveform "x".
% ( see IDIAP RR 06-09 ). This applies CHN, then USS,
% and returns the filtered waveform "chnuss_x".
%
% Here, CHN is applied in the frequency domain.
% In "chnuss_signal_main2.m", CHN is applied in the time domain.
%
% Mandatory parameters:
%   par.block_size_frames: size of block, in frames (e.g. 100).
%   par.do_chn = 0 or 1
%   par.do_uss = 0 or 1
%   par.preemp = 0.97

% Default value for par.chn_alpha: 0.2
% 
% You can set some other parameters, to see them do:
%   dbtype chnuss_signal_main2.
%
%
% Guillaume Lathoud - lathoud@idiap.ch - 2006


% Check parameters
  
  if nargin < 2
    error( [ mfilename ' needs at least 2 input arguments.' ] );
  end

  check_param( { 'block_size_frames', 'do_chn', 'do_uss', 'preemp' }, fieldnames( par ) );
  
  par_default = [];

  par_default.framelength = 512 ;    % WINDOWSIZE converted to samples instead of 100 nanosecond crap
  par_default.frameshift = [] ;      % TARGETRATE (ditto from WINDOWSIZE)
  par_default.fs = 16000 ;            % sampling frequency in Hz

  par_default.past_nframes_max   = par.block_size_frames;
  par_default.future_nframes_max = par.block_size_frames;
  
  par_default.chn_alpha = 0.2;   % For the channel estimation (IDIAP RR 06-09)
  
  par_default.sanity_check = 1;

  par_default.debug = 0; % if 1 then show some figures and make some pauses
  
  par_default.verbose = 1;

  par = fill_default( par, par_default );
  
  if isempty( par.frameshift )
    par.frameshift = round( par.framelength / 2 );
  end
  
  if par.verbose
    disp( '--------------------------------------------------' );
    disp( [ mfilename ' par:' ] );
    disp( par );
  end

  if par.debug
    h = [ figure figure figure figure ];
  end
  
  % Define Hamming window
  ham_win = diag( sparse( hamming( par.framelength, 'periodic' ) ) );
  
  % initialization
  
  x = x(:);
  
  length_y = length( x );

  start_sample = 1;
  stop_sample  = length_y;

  num_frames = 1 + floor( ( length_y - par.framelength ) / par.frameshift ) ;
  
  % Where we will store the output
  chnuss_x      = zeros( length_y, 1 );
  
  par.block_size_frames = min( num_frames, par.block_size_frames );
  
  if par.verbose
    chrono = chrono_start( num_frames );
  end
  
  % Initialization: Online channel normalization
  chn_obj.past_data = [];

  % Initialization: Online Unsupervised Spectral Subtraction
  if par.do_uss
    rse_uss( 'init', [] );
  end

  
  % frame by frame processing

  start_frame = 1;
  stop_frame  = num_frames;

  for i = start_frame:par.block_size_frames:stop_frame
    
    % ( 1 ) Extract spectrum for this frame, if not already done
    %       (this is done per blocks of frames)
    %       This allows the user to select offline / blockwise online / frame-by-frame online 
    %       processing, by simply changing "par.block_size_frames".
    
    a_block.start_frame  = i;
    a_block.stop_frame    = min( num_frames, a_block.start_frame + par.block_size_frames - 1 );
    a_block.block_size_frames      = 1 + a_block.stop_frame - a_block.start_frame;
    
    % "0" because we need an extra previous sample for pre-emphasis
    a_block.start_sample =                0 + ( a_block.start_frame - 1 ) * par.frameshift;
    
    a_block.stop_sample   = par.framelength + ( a_block.stop_frame - 1 ) * par.frameshift;
    
    % Extract the frames of this block
    % "max( 1," because the extra previous sample is not available at the beginning
    % -> the first sample of the first block is duplicated.
    a_block.y            = x( max( 1, a_block.start_sample:a_block.stop_sample ) );
    
    a_block.frame_mat    = buffer( a_block.y, 1 + par.framelength, 1 + par.framelength - par.frameshift, 'nodelay' );
    
    if par.sanity_check
      % Sanity check
      if ~all( size( a_block.frame_mat ) == [ 1 + par.framelength  a_block.block_size_frames ] )
	error( [ mfilename ' is insane!' ] );
      end
    end
    
    % zero mean
    % ZMEANSOURCE = T
    a_block.frame_mat = a_block.frame_mat - repmat( mean( a_block.frame_mat, 1 ), size( a_block.frame_mat, 1 ), 1 );
    
    % preemphasise NOT according to formula in book
    % -> let's not introduce unnecessary noise in the signal
    % PREEMCOEF = P
    a_block.frame_mat = a_block.frame_mat( 2:end,: ) - par.preemp * a_block.frame_mat( 1:end-1,: );
    
    % window using Hamming window
    % USEHAMMING = T
    a_block.frame_mat = ham_win * a_block.frame_mat;
    
    % do zero-padded fft
    a_block.f_cplx_frame_mat = fft( a_block.frame_mat , 2 * par.framelength, 1 );
    a_block.f_frame_mat      = abs( a_block.f_cplx_frame_mat );

    if par.debug
      figure( h( 1 ) ); clf; imagesc( log( a_block.f_frame_mat ) );
      title( strrep( [ mfilename ': log mag spectrum before filtering (one block).' ], '_', '\_' ) );
    end
    
    % ( 2.1 ) Channel normalization
    
    % Get the "future" data for channel normalization
    
    tmp   = a_block.f_frame_mat;
    
    tmp( find( tmp == 0 ) ) = +Inf;
    
    tmp   = sort( tmp, 2 );
    
    % 20 % lower non-zero magnitudes
    % This is directly inspired from Prof. Rainer Martin's work,
    % although he was using it for noise level estimation,
    % while we are using it for channel estimation.
    
    n     = size( tmp, 2 );
    tmp   = tmp( :, 1:1 + round( ( n - 1 ) * par.chn_alpha ) );
    
    % Check past and future have about the same sizes
    % and are not too large
    
    n     = size( tmp, 2 );
    n_max = par.future_nframes_max; 
    
    if n > n_max
      tmp = tmp( :, 1 + round( ( n - 1 ) * ( 0.5:n_max ) / n_max ) );
    end
    
    % Make sure past is not larger than future
    
    n     = size( chn_obj.past_data, 2 );
    n_max = par.past_nframes_max;

    if n > n_max
      chn_obj.past_data = sort( chn_obj.past_data, 2 );
      chn_obj.past_data = chn_obj.past_data( :, 1+round( ( n-1 ) * ( 0.5:n_max ) / n_max ) );
    end

    % Data on which we potentially estimate the channel

    tmp = [ chn_obj.past_data  tmp ];
    
    if ~par.do_chn

      % No channel normalization

      channel = ones( size( tmp, 1 ), 1 );

    else
      
      % Estimate the channel using "future" and "past" data
      
      channel = NaN + ones( size( tmp, 1 ), 1 );
      
      for qq = 1:( 1 + par.framelength )
	
	channel( qq ) = exp( mean( log( tmp( qq, find( ( ~isnan( tmp( qq, : ) ) ) & ( ~ismember( tmp( qq,: ), [ 0 +Inf -Inf ] ) ) ) ) ) ) );
	
      end
      
      % Smooth the channel along the frequency axis
      
      ind = 1:1+par.framelength;
      channel( ind )   = exp( running_mean( log( channel( ind ) ), 5 ) );
      
      % Replicate to negative freqs
      
      channel( par.framelength+2:end ) = channel( par.framelength:-1:2 );

      % In some special case we don't have enough data
      % (e.g. all-zero signal due to a telephone VAD)
      
      if all( isnan( channel ) )
	channel = ones( size( channel ) );
      end
      
      % Some other cases
      
      if any( isnan( channel ) )
	channel( find( isnan( channel ) ) ) = exp( mean( log( channel( find( ~isnan( channel ) ) )  )));
      end
      
      % Sanity check
      
      if any( isnan( channel ) )
	error( [ mfilename ': channel is NaN-insane.' ] );
      end
      
      if par.debug
	figure( h( 2 ) ); clf; plot( log( channel ) ); title( strrep( [ mfilename ': log channel (one block).' ], '_', '\_' ) );
      end
      
    end

    % Update the "past" data,
    % which will be used while processing the next block
    
    chn_obj.past_data = sort( tmp, 2 );
    
    n = size( chn_obj.past_data, 2 );
    
    n_max = par.past_nframes_max;
    
    if n > n_max
      
      chn_obj.past_data = chn_obj.past_data( :, 1 + round( ( n - 1 ) * ( 0.5:n_max ) / n_max ) );
      
    end
    

    % ( 2.2 ) Unsupervised Spectral Subtraction (USS)
    % on the channel-normalized magnitude spectrum

    a_block.chn_f_frame_mat = diag( sparse( 1 ./ ( eps + channel ) ) ) * a_block.f_frame_mat;

    if par.debug
      figure( h( 3 ) ); clf; imagesc( log( abs( a_block.chn_f_frame_mat ) ) );
      title( strrep( [ mfilename ': log mag spectrum after chn (one block).' ], '_', '\_' ) );
    end

    if ~par.do_uss

      % USS not activated

      a_block.chnuss_f_frame_mat = a_block.chn_f_frame_mat;

    else

      % USS activated

      a_block.chnuss_f_frame_mat = rse_uss( a_block.chn_f_frame_mat, par );

    end
    
    if par.debug
      figure( h( 4 ) ); clf; imagesc( log( abs( a_block.chnuss_f_frame_mat ) ) );
      title( strrep( [ mfilename ': log mag spectrum after chn and uss (one block).' ], '_', '\_' ) );
    end
    
    % ( 3 ) Filter that portion of the signal
    
    filter_mat = a_block.chnuss_f_frame_mat ./ ( eps + a_block.f_frame_mat );
    
    % Filter the original signal, using overlap-add
    
    a_block.chnuss_x = overlap_add( a_block.f_cplx_frame_mat, filter_mat, par.framelength, par.frameshift, par.fs, par.preemp );

    % Add a smooth transition at the beginning and end of each block.
    % We use the fact that cos^2+sin^2 = 1

    if a_block.start_frame > start_frame
      
      % End of previous block: smooth decay

      ind_1 = a_block.start_sample + 1;   % +1 because of pre-emphasis (see above)
      ind_2 = a_block.previous_stop_sample;

      if ind_2 > ind_1
	
	kk = ind_1:ind_2;
	kk = kk(:);
	
	chnuss_x( kk ) = chnuss_x( kk ) .* ( cos( pi / 2 * ( kk - ind_1 ) / ( ind_2 - ind_1 ) ) .^ 2 );
	
      end

      % Beginning of current block: smooth increase
      
      ind_1 = 1;
      ind_2 = 1 + a_block.previous_stop_sample - a_block.start_sample - 1;  % -1 because of pre-emphasis (see above)
      
      if ind_2 > ind_1

	kk = ind_1:ind_2;
	kk = kk(:);
	
	a_block.chnuss_x( kk ) = a_block.chnuss_x( kk ) .* ( sin( pi / 2 * ( kk - ind_1 ) / ( ind_2 - ind_1 ) ) .^ 2 );
    	
      end

    end
    
    a_block.previous_stop_sample = a_block.stop_sample;
    
    % Add the result
    
    ind = a_block.start_sample+1:a_block.stop_sample;
    chnuss_x( ind ) = chnuss_x( ind ) + a_block.chnuss_x;
    
    if par.verbose
      chrono = chrono_check( chrono, i );
    end
    
    if par.debug
      disp( [ mfilename ': debug: press return to continue.' ] );
      pause;
    end

  end  % Loop through blocks
  
  if any( isnan( chnuss_x(:) ) )
    error( [ mfilename ' is NaN-insane.' ] );
  end
  
  if any( isinf( chnuss_x(:) ) )
    error( [ mfilename ' is Inf-insane.' ] );
  end

  if par.verbose
    
    chrono_stop( chrono );
    disp( [ mfilename ': done' ] );
    
    disp( ' ' );
    
  end
