function mag = chn( mag, par )
  
% Implementation of CHannel Normalization (CHN)
%
% See IDIAP RR 06-09.
%
% By Guillaume Lathoud 2006
% lathoud@idiap.ch

  if nargin < 2
    error( [ mfilename ' needs 2 input arguments.' ] );
  end
  
  % Where we'll store the "forgetting" memory
  persistent past_data;

  if ischar( mag )
    
    % Initialization
    past_data = [];
    return;
    
  end

  % Check input parameters
  required_fields = { 'future_nframes_max', 'past_nframes_max', 'chn_prop' };
  check_param( required_fields, fieldnames( par ) );
  
  % Make sure this is not a complex spectrum
  if ~isreal( mag )
    
    error( [ mfilename ' needs an array of (real) magnitude values for "mag".' ] );
    
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 1 ) Update the data representation
  %       (similar to a forgetting factor)
  
  % ( 1.1 ) Reduce the size of the "past" data
  
  past_data = sort( past_data, 2 );
  
  n     = size( past_data, 2 );
  n_max = par.past_nframes_max;

  if n > n_max

    past_data = past_data( :, round( n * ( 0.5:n_max ) / n_max ) );

  end
  

  % ( 1.2 ) Reduce the size of the "future" data
  
  sorted_mag = sort( mag, 2 );
  
  n     = size( sorted_mag, 2 );
  n_max = par.future_nframes_max;

  if n > n_max
    
    sorted_mag = sorted_mag( :, round( n * ( 0.5:n_max ) / n_max ) );

  end
  
  % ( 1.3 ) Update
  
  past_data = [ past_data   sorted_mag ];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 2 ) CHannel Normalization (CHN)
  
  % ( 2.1 ) Within each frequency band, the channel is estimated 
  % as the geometric mean of the lowest "par.chn_prop" magnitude values.

  channel = repmat( NaN, size( mag, 1 ), 1 );
  
  for a = 1:numel( channel )

    ind = find( ~isnan( past_data( a,: ) ) & ~isinf( past_data( a,: ) ) & ( past_data( a,: ) > 0 ) );
    
    if ~isempty( ind )
      
      % Pick the lowest magnitude values (most likely to be noise, see Rainer Martin Eurospeech 93)
      tmp = sort( past_data( a, ind ) );
      tmp = tmp( 1:max( 1, round( numel( tmp )  *  par.chn_prop ) ) );

      % Compute the geometric mean
      channel( a ) = exp( mean( log( tmp ) ) );

    end
    
  end

  % Smooth out outliers with a 5-point running mean
  channel = running_mean( channel, 5 );

  % Deal with very rare case
  if any( isnan( channel ) )
    channel( find( isnan( channel ) ) ) = exp( mean( log( [ 1; channel( find( ~isnan( channel ) ) ) ] ) ) );
  end
  
  if par.debug
    figure; semilogy( channel );
    title( [ mfilename ' - magnitude channel' ] );
    keyboard
  end

  % ( 2.2 ) CHN
  
  mag = full( diag( sparse( 1 ./ ( eps + channel ) ) ) * mag );
  