function mag = rse_uss( mag, par )
  
% Implementation of fitting of Rayleigh + Shifted Erlang (RSE)
% and Unsupervised Spectral Subtraction (USS)
% (ASRU 2005 - Lathoud, Magimai.-Doss and Bourlard)
% (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', 'rse_band', 'rse_em_par' };
  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, using only data within a specified band
  % ( use rse_band = [ -Inf +Inf ] for fullband )
  
  band_ind = max( 1, par.rse_band( 1 ) ):min( size( sorted_mag, 1 ), par.rse_band( 2 ) );
  
  past_data = [ past_data   sorted_mag( band_ind, : ) ];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 2 ) RSE-USS
  
  % ( 2.1 ) The noise level is estimated by fitting a Rayleigh + Shifted Erlang model,
  % irrespective of time or frequency (RSE)
  
  ind = find( ~isnan( past_data(:) ) & ~isinf( past_data(:) ) & ( past_data(:) > 0 ) );
  
  % By default, result of the filtering will be silence (hence the default infinite noise level)
  noise_level = +Inf;
  
  if ~isempty( ind )

    % Usual case: there is at least one proper magnitude value

    % Reduce the data, we don't need many samples
    tmp = sort( past_data( ind ) );
    n = numel( tmp );
    n_max = 100;
    if n >= n_max
      tmp = tmp( round( n * ( 0.5:n_max ) / n_max ) );
    end
    
    % Fit the RSE model
    raylsherl = fit_raylsherl( tmp, par.rse_em_par );
    
    if isempty( raylsherl )
      
      % Failure
      noise_level = +Inf;
      
    else
      
      % Success
      % Get the estimated noise level
      noise_level = raylsherl.erlang_t;
      
    end
    
  end

  
  % ( 2.2 ) Unsupervised Spectral Subtraction (USS)
  
  mag = max( 1, mag / noise_level );
  