function  a_struct = FASTTDE_detect_locate( a_struct, ssm, PAR )

% function  a_struct = FASTTDE_detect_locate( a_struct, ssm, PAR )
%
% Process a multichannel recording in an online fashion,
% by iterating through the following steps:
% A. instantaneous sector-based multisource detection (IDIAP RR 05-52).
% B. instantaneous multisource localization (within each active sector, based on Knapp 76 and Brandstein 95).
% C. short-term clustering of instantaneous location estimates (IDIAP RR 04-14).
% D. optionally, extraction of sector-based MFCCs (for speech/non-speech, IDIAP RR 06-26).
% E. optionally, extraction of single-channel MFCCs (for speaker clustering, IDIAP RR 06-26).
%
% "FASTTDE" means that we use a direct localization method, based on GCCPHAT
% Time-Delay Estimation.
%
% Moreover, a fast implementation is given for sector-based multisource detection.
% It computes SAM-SPARSE-MEAN values (IDIAP RR 04-67) from mic. array signals.
% "get_ssm_activeness_c.c" is required.                                       
% Use "mex get_ssm_activeness_c.c" to compile it on your machine. 
% The -ffast-math option does not seem to be a problem with it.                                
%   A "2 Diracs + 2 Gammas + 1 Shifted Erlang" model is then used,
% as in IDIAP RR 05-52. Shifted Erlang is more stable than Shifted Rice,
% for limited amounts of data.
%   The fitted model is used with two different unsupervised thresholds
% (PAR.DETECTION.FAR_T and PAR.DETECTION.FRR_T), and the results are
% then combined using PAR.DETECTION.SPACE_TIME_DILATION. The idea is
% that a segment should have at least one frame above the conservative
% threshold (FAR_T).
% 
% The localization within each active sector (step B.) is made through 
% GCC-PHAT Time Delay Estimation (FASTTDE). For each square subarray,
% a direction is estimated using Branstein's thesis (1995, Section 7.2).
% It relies on an upsampled time-domain GCC-PHAT.
%
% Short-term clustering (C.) is applied in azimuth space only.
% (elevation and radius are not reliable with a UCA).
% (MATLAB implementation).
%
% if PAR.DO_SBMFCC == 1:
% Sector-based MFCCs (D.) are extracted from a sector-based spectrum,
% itself obtained by zeroing magnitude at frequencies where a sector
% is not dominant. "Dominance" is determined as in IDIAP RR 04-67.
% To avoid unnecessary dynamics in the cepstrum domain, we then floor
% the spectrum as in Unsupervised Spectral Subtraction (IDIAP RR 06-09).
%
% if PAR.DO_MFCC == 1:
% Single-channel MFCCs (E.) are basic MFCCs extracted from the raw,
% unmodified signal of channel PAR.MFCC_CHANNEL.
%
% Note: useful things related to MEX and C compilation
% (1) "mex -setup" to create a configuration file 
%     (e.g. "~/.matlab/R13/mexopts.sh")
% (2) In that file, you can set useful options such as
%                COPTIMFLAGS='-O3 -DNDEBUG -march=pentium4 -ffast-math'
%
% By Guillaume Lathoud, May 2006
% lathoud@idiap.ch
%
% Thanks to Bertrand Mesot for advices on MEX.
  
  if nargin < 3
    error( [ mfilename ' needs 3 input arguments.' ] );
  end
  

  % To pass large pieces of data efficiently,
  % to and from the various subfonctions
  % (otherwise their copied by MATLAB)
  global ws  % Stands for "workspace"
  ws = [];
  
  % To communicate with the short-term clustering
  % "get_stc_past_future" and "update_stc_past_future"
  global sphloc_mat
  sphloc_mat = [];

  % Paranoid setting
  dbstop error
  dbstop warning

  % To count time
  start_clock = clock;
  
  a_struct.total_clock     = 0;  % Total time spent in this function.
  a_struct.total_fasttde_clock = 0;  % Total time spent doing GCC-PHAT TDE-based localization
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 0 ) Check parameters
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % Mandatory
  
  check_mandatory( a_struct, ssm, PAR );
  
  % Optional: fill with default values

  PAR = check_optional( PAR, ssm );

  % Details: framelength, shift etc.

  a_struct = get_details( a_struct, ssm, PAR );


  if PAR.VERBOSE
    disp( ' ' );
    disp( [ mfilename ' parameters PAR:' ] );
    disp( PAR );
    disp( [ mfilename ' parameters PAR.DETECTION:' ] );
    disp( PAR.DETECTION );
    disp( [ mfilename ' parameters PAR.STC:' ] );
    disp( PAR.STC );
    if PAR.DO_MFCC
      disp( [ mfilename ' parameters PAR.USS:' ] );
      disp( PAR.USS );
      disp( [ mfilename ' parameters PAR.MFCC:' ] );
      disp( PAR.MFCC );
    end

    disp( sprintf( 'framelength: %d samples,  frameshift: %d samples,  fs: %d Hz,  c: %6.2f m/s,  fft_framelength: %d samples', a_struct.framelength, a_struct.frameshift, ssm.fs, ssm.c, a_struct.fft_framelength ) );
    disp( sprintf( 'start_frame: %d,  stop_frame: %d,  block size: %d frames,  init block size: %d frames', a_struct.start_frame, a_struct.stop_frame, a_struct.block_size_frames, a_struct.init_block_size_frames ) );
  end
  
  % For compatibility with the C implementation "get_ssm_activeness_c.c"

  ssm.pairs = uint16( ssm.pairs );
  
  % For easy access from within subfunctions

  ws.a_struct = a_struct;
  ws.ssm = ssm;
  ws.PAR = PAR;
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 1 ) Initialize the detection model
  % using the beginning of the recording.

  if PAR.VERBOSE
    disp( ' ' );
    disp( sprintf( [ mfilename ': initializing the detection model with %d frames...' ], a_struct.init_block_size_frames ) );
  end

  % ( 1.1 ) Get the frequency-domain signals 
  % They are stored into "ws.block.X".
  %
  % Also, GCC-PHAT is computed and stored into
  % "ws.block.ssm_observed_gccphat".
  
  if PAR.DETECTION.DEBUG
    disp( [ '[          Extracting GCC-PHAT from waveforms...          ]' ] );
  end
  
  ws.block.start_frame = a_struct.start_frame;
  ws.block.stop_frame  = min( ws.block.start_frame + a_struct.init_block_size_frames - 1, a_struct.stop_frame );
  
  get_observed_gccphat;
  
  % ( 1.2 ) Get SSM activeness from complex domain GCC-PHAT
  % measurements store it into "ws.block.ssm_activeness"
  % ( as well as the underlying binary decision
  % into "ws.block.ssm_binary_decision" ).
  %
  % SSM stands for SAM-SPARSE-MEAN (IDIAP RR 04-67).

  if PAR.DETECTION.DEBUG
    disp( [ '[          Extracting SSM activeness from GCC-PHAT...          ]' ] );
  end
  
  get_ssm_activeness;
  
  % ( 1.3 ) Reduce the data to a fixed size
  % This creates a field "ws.block.reduced_ssm_activeness".
  
  if PAR.DETECTION.DEBUG
    disp( [ '[          Reducing data size...          ]' ] );
  end
  
  ws.block.previous_ssm_activeness = [];

  reduce_ssm_data;
  
  % ( 1.4 ) Unsupervised fit of the detection model
  % on "ws.block.reduced_ssm_activeness"
  %
  % (IDIAP RR 05-52)
  
  if PAR.DETECTION.DEBUG
    disp( [ '[          Fitting a model...          ]' ] );
  end
  
  fit_model;
  
  % ( 1.5 ) Unsupervised threshold selection,
  % thresholds stored in "ws.block.far_threshold" and "ws.block.frr_threshold".
  %
  % (IDIAP RR 05-52)

  if PAR.DETECTION.DEBUG
    disp( [ '[          Selecting thresholds...          ]' ] );
  end
  
  select_thresholds;
  

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 2 ) Initialize the "same/different" local models
  % for short-term clustering (IDIAP RR 04-67).
  %
  % Construct all possible graphs with up to PAR.STC.N_FUTURE_SAMPLES vertices.
  
  if PAR.DO_STC
    
    if PAR.VERBOSE
      disp( ' ' );
      disp( [ mfilename ': fixed initialization of the short-term clustering local models...' ] );
    end
    
    ws.stc = [];

    % ( 2.1 ) Local models
    
    % Our model uses zero-centered Gaussians,
    % for azimuth (az).
    %
    % We initialize "as if" the variances did not depend on the time interval
    %
    % The model is described by two matrices of STD DEV values.
    % ( zero-mean is assumed ):
    % - "same":    Gaussian with 66% contained between -max, +max.
    % - "notsure": uniform on ( [ -pi, -max [ union ] +max, +pi ] ).
    %
    % Both "same" and "notsure" have only one row, describing azimuth space:
    % - column 1: delta_t = 0
    % - column 2: delta_t = 1
    % ...
    % - column 1 + PAR.STC.T_SHORT:  delta_t = PAR.STC.T_SHORT
    
    a = ones( 1, PAR.STC.T_SHORT + 1 );
    
    m = PAR.STC.MAX_DAZ;

    ws.stc.model.notsure = [ a *  sqrt( ( pi^2 + pi * m + m^2 ) / 3 ) ];

    ws.stc.model.same    = [ a *  m / 2 ];
    
    % Equal priors
    
    ws.stc.model.p_notsure = 0.5 * a;
    ws.stc.model.p_same    = 1 - ws.stc.model.p_notsure;
    
    % start without data
    
    ws.block.reduced_stc_data = {};

    % Mark the last location estimate that has a STC tag
    % This is a column index in sphloc_mat
    ws.stc.last_loc_est = [];  % None so far
    ws.stc.last_frame  = ws.a_struct.start_frame - 1; % No frame so far

    % Mark the new samples/frames available
    ws.stc.n_new_sample  = 0;  % None so far
    ws.stc.new_frame     = ws.a_struct.start_frame;

    % Also initialize the tag counter.
    % We'll increment it each time a cluster is created.
    ws.stc.last_tag      = 0;

    % ( 2.2 ) Construct all possible graphs with up to PAR.STC.N_FUTURE_SAMPLES vertices.
    
    ws.stc.all_graphs = [];
    
    for a = 1:PAR.STC.N_FUTURE_SAMPLES
      
      if a < 2
	ws.stc.all_graphs( a ).graph_list = construct_graph_list_pruned( a, [], [], [], PAR.STC.N_GRAPHS_MAX );
      else
	ws.stc.all_graphs( a ).graph_list = construct_graph_list_pruned( a, ws.stc.all_graphs( a-1 ).graph_list, [], [], PAR.STC.N_GRAPHS_MAX );
      end
      
    end
    
    if PAR.STC.VERBOSE
      disp( ' ' );
      disp( sprintf( [ mfilename ': constructed a total number of %d possible partitions of up to %d elements' ], length( ws.stc.all_graphs( end ).graph_list ), PAR.STC.N_FUTURE_SAMPLES ) );
    end
    
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 3 ) Process the multichannel recording
  % in an online manner. See the main comment above.
  %
  % Everything is online, frame-by-frame:
  % to speedup things, we access the files block per block,
  % but EACH frame is processed individually.
  % -> Model parameters are updated at the end of each block.
  
  % Where we'll store the detection decision
  % ( zeroes and ones )
  % The column index is the absolute frame index from the start of the recording.
  
  ws.ssm_detection_decision = logical( sparse( a_struct.n_sectors, a_struct.n_frames ) );
  
  % Where we will store the matrix containing location estimates.
  %
  % Format of this matrix: each column contains one location estimate.
  % Note that there can be zero, one or more location estimate per time frame.
  % Therefore, the columns in "sphloc_mat" have no relation with time frames.
  %   row 1: absolute frame index (integer), i.e. from the start of the recording.
  %   row 2: sector index (integer).
  %   row 3: azimuth in radians (very reliable with a Uniform Circular Array or UCA).
  %   row 4: elevation in radians (not very reliable with a UCA).
  %   row 5: radius in meters (not reliable at all with a UCA).
  %   row 6: log posterior probability of activity (following IDIAP RR 05-52).
  %   row 7: (optional) tag given by short-term clustering (NaN if none).

  sphloc_mat = zeros( 7, 0 );
  
  if PAR.DO_SBMFCC
    
    ws.PAR.MFCC_DIM =  ws.PAR.MFCC.N + ws.PAR.MFCC.DO_C0;

    % Where we will store the matrix
    % containing sector-based MFCC's of each active (time frame, sector)
    % ( conservative MFCC's -> loss of information )
    
    %   row 1: frame index (integer)
    %   row 2: sector index (integer)
    %   row 3 to end: MFCC vector
    global sbmfcc_mat
    
    sbmfcc_mat = zeros( 2 + ws.PAR.MFCC_DIM, 0 );
    
  end

  if PAR.DO_MFCC

    ws.PAR.MFCC_DIM =  ws.PAR.MFCC.N + ws.PAR.MFCC.DO_C0;

    % Where we'll store MFCCs extracted from the raw signal
    % of one microphone channel.
    
    global mfcc_mat
    
    mfcc_mat = zeros( 1 + ws.PAR.MFCC_DIM, 0 );
    
  end


  % Prepare a exchange structure for one_array_fasttde_dirloc_c.c
  % (MEX implementation of TDE + direction estimation)

  global fasttde_wsp

  fasttde_wsp = [];
  fasttde_wsp.i_i_pair_list             = uint16( full( ws.a_struct.fasttde.i_i_pair_list ).' );  % transpose
  
  for a = 1:ws.a_struct.fasttde.nsubarray
    fasttde_wsp.square_subarray_Pmat_list{ a } = full( ws.a_struct.fasttde.square_subarray_Pmat_list{ a } );
  end
  
  fasttde_wsp.up_rowzero                = uint16( full( ws.a_struct.fasttde.up_rowzero ) );
  fasttde_wsp.upfactor                  = uint16( full( ws.a_struct.fasttde.upfactor ) );
  fasttde_wsp.maxdelay                  = full( ws.a_struct.fasttde.maxdelay );

  %
  
  if PAR.VERBOSE
    disp( ' ' );
    disp( [ mfilename ': processing frames in an online fashion...' ] );
    chrono = chrono_start( a_struct.n_frames, 30 );
  end
  
  % To measure actual time spent in the processing.
  % This is meaningless if you are running other processes at the same time.
  % 
  % Global time:
  a_struct.init_clock = etime( clock, start_clock );
  %
  % Time for each step: for time measurement,
  % we have splitted each iteration in 12 different steps.
  a_struct.total_step_clock = zeros( 1, 12 );
  
  for i_frame = 1:a_struct.n_frames
    
    % ---------- CLOCK STEP LIMIT
    tmp_clock = clock;

    % To make data reading faster, 
    % we do as much processing as possible per block of frames,
    % without violating the "frame-by-frame", "online" constraint
    % mentioned above.
    
    is_new_block = ~mod( i_frame - 1, a_struct.block_size_frames );
    
    % Absolute frame index from the beginning of the file (starting with 1)
    
    frame_ind    = a_struct.start_frame + i_frame - 1;

    % ( 3.1 ) Get the frequency-domain signals
    % They are stored into "ws.block.X".
    %
    % Also, frequency-domain GCC-PHAT is computed and stored into
    % "ws.block.ssm_observed_gccphat".
    
    if is_new_block
      
      ws.block.start_frame = frame_ind;
      ws.block.stop_frame  = min( ws.block.start_frame + a_struct.block_size_frames - 1, a_struct.stop_frame );
      
      if PAR.DETECTION.DEBUG
	disp( sprintf( [ '[          Extracting GCC-PHAT from waveforms (frames %d to %d)...          ]' ], ws.block.start_frame, ws.block.stop_frame ) );
      end
      
      get_observed_gccphat;
 
    end

    % Relative frame index from the beginning of the current block
    
    block_frame_ind = 1 + frame_ind - ws.block.start_frame;
    

    a_struct.total_step_clock(1) = a_struct.total_step_clock(1) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT
    tmp_clock = clock;

    % ( 3.2 ) Extract the SSM activeness
    
    if is_new_block
      
      if PAR.DETECTION.DEBUG
	disp( [ '[          Extracting SSM activeness from GCC-PHAT...          ]' ] );
      end
      
      get_ssm_activeness;
            
    end
    
    a_struct.total_step_clock(2) = a_struct.total_step_clock(2) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;
    
    % ( 3.3 ) Apply the detection model
    
    if is_new_block
      
      if PAR.DETECTION.DEBUG
	disp( [ '[          Computing posteriors...          ]' ] );
      end

      [ tmp, sector_log_post ] = compute_post_2d2g1se( ws.block.model, ws.block.ssm_activeness );
  
      % We are only interested in the sector-level posterior of activity
      % (IDIAP RR 05-52)
      
      ws.block.ssm_log_post_act = sector_log_post( :,:,3 );
      
      % Take the detection decision by comparing the posterior
      % to thresholds, store the decision.

      if PAR.DETECTION.DEBUG
	disp( [ '[          Taking the detection decision...          ]' ] );
      end

      take_detection_decision;
      
      % If needed, compute time-domain GCC-PHAT
      % We'll use it for time-delay estimation
      
      tmp = clock;  % This is part of time-delay estimation

      compute_timedomain_gccphat;
      
      a_struct.total_fasttde_clock = a_struct.total_fasttde_clock + etime( clock,  tmp );
      
    end
    
    a_struct.total_step_clock(3) = a_struct.total_step_clock(3)  + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    
    % ( 3.4 ) 3-D localization within each active sector
    % Relies on an MATLAB implementation of GCC-PHAT time-delay estimation
    
    active_sector_list = find( ws.ssm_detection_decision( :, frame_ind ) );

    % Limit the search
    
    if numel( active_sector_list ) > PAR.DETECTION.N_SECTORS_MAX
      
      [ tmp, ind ] = sort( -ws.block.ssm_log_post_act( active_sector_list, block_frame_ind ) );
      
      active_sector_list = active_sector_list( ind( 1:PAR.DETECTION.N_SECTORS_MAX ) );
      
    end
    
    % If necessary, fetch the time-domain GCCPHAT
    
    if numel( active_sector_list ) > 0

      fasttde_wsp.up_gccphat = full( ws.block.up_gccphat( :, :, block_frame_ind ) );

    end

    % Locate within each active sector

    a_result = [];
    
    for sector_id = active_sector_list(:).'
      
      if PAR.DETECTION.DEBUG
	disp( sprintf( [ mfilename ': frame %d, localization within sector %d...' ], frame_ind, sector_id ) );
      end
      
      % Do the GCCPHAT TDE-based localization within this sector.
      % It may fail, in which case [] are returned for both "sub_az" and "sub_al"
      
      tmp = clock;
      
      if 1
	
	% Call the C implementation

	% Set the inputs
	fasttde_wsp.az_min         = full( ws.ssm.sector( sector_id ).az_min );
	fasttde_wsp.az_max         = full( ws.ssm.sector( sector_id ).az_max );
	fasttde_wsp.el_min         = full( ws.ssm.sector( sector_id ).el_min );
	fasttde_wsp.el_max         = full( ws.ssm.sector( sector_id ).el_max );
	fasttde_wsp.up_index_min   = uint16( full( ws.a_struct.fasttde.up_index_min( :, sector_id ) ) );
	fasttde_wsp.up_index_max   = uint16( full( ws.a_struct.fasttde.up_index_max( :, sector_id ) ) );
	
	% Allocate memory for the outputs
	fasttde_wsp.az = 0;
	fasttde_wsp.el = 0;
	fasttde_wsp.success = 0;

	one_array_fasttde_dirloc_c;
	
	if fasttde_wsp.success

	  az = fasttde_wsp.az;
	  el = fasttde_wsp.el;

	else

	  az = [];
	  el = [];

	end

      else
	
	% MATLAB implementation, kept here to explain what the C implementation does
	[ az, el ] = one_array_fasttde_dirloc( block_frame_ind, sector_id );
	
      end
      
      a_struct.total_fasttde_clock = a_struct.total_fasttde_clock + etime( clock,  tmp );
      
      % Success?

      if ( ~isempty( az ) ) & ( ~isempty( el ) )
	
	% Yes
	
	v        = [ frame_ind;
		     sector_id;
		     az;
		     el ];
	
	% Radius unknown
	v( end + 1, : ) = NaN;
	
	% Also store the log posterior probability of activity
	v( end+1,: ) = ws.block.ssm_log_post_act( sector_id, block_frame_ind );
	
	% Unknown cluster tag at this point
	v( end+1,: ) = NaN;
	
	if ws.PAR.DETECTION.DEBUG
	  disp( [ '[ ' sprintf( ' (frame_ind:%d  sector_id:%d  az:%.3f  el:%.3f  r:%.3f  log_post_act:%g  cluster_tag:%d) ', v ) ' ]'  ] );
	end
	
	% Stack it
	a_result = [ a_result  v ];
	
      end
      
    end  % of loop through sectors

    n_points = 0;
    
    if ~isempty( a_result )
    
      % Store all location estimates for this frame
      sphloc_mat = [ sphloc_mat   a_result ];
      
      n_points   = size( a_result, 2 );

    end

    a_struct.total_step_clock(4) = a_struct.total_step_clock(4)  + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    % ( 3.5 ) Short-term clustering of locations

    if PAR.DO_STC
      
      % Regular treatment: update the counts
      %
      % "new sample" means "does not have a cluster tag"

      ws.stc.n_new_sample  = ws.stc.n_new_sample + n_points;
      ws.stc.new_frame     = frame_ind + 1;

      while ws.stc.n_new_sample >= ws.PAR.STC.N_FUTURE_SAMPLES
	
	% If too much time has passed and/or if enough data has accumulated,
	% then do short-term clustering
	
	do_stc = ( ws.stc.n_new_sample >= ws.PAR.STC.N_FUTURE_SAMPLES );
	do_stc = do_stc | ( ws.stc.new_frame > ws.stc.last_frame + ws.PAR.STC.T_SHORT );
	
	% Also do it when reaching the last frame
	
	do_stc = do_stc | ( i_frame == ws.a_struct.n_frames );
	
	% ... but in all cases, we need at least one new location estimate!
	
	do_stc = do_stc   &   ( ws.stc.n_new_sample > 0 );
	
	
	if do_stc
	  
	  if PAR.STC.VERBOSE
	    disp( [ '[          Short-term clustering of location estimates...          ]' ] );
	  end
	  
	  % Get the short-term past and future
	  % This modifies: "ws.stc.past" and "ws.stc.future"
	  get_stc_past_future;
	  
	  % Apply short-term clustering
	  % This modifies the tags in "ws.stc.future"
	  short_term_clustering;
	  
	  % Store the result
	  % This modifies: "sphloc_mat", "ws.stc.last_loc_est"
	  update_stc_past_future;
	  
	else
	  
	  % Nothing left to do at the moment,
	  % move on to the next frame.

	  break;
	  
	end
	
      end
      
    end
    
    a_struct.total_step_clock(5) = a_struct.total_step_clock(5) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    % ( 3.6 ) Reduce the sector activity data to a fixed size
    
    if frame_ind == ws.block.stop_frame
      
      if PAR.DETECTION.DEBUG
	disp( [ '[          Reducing sector activity data size...          ]' ] );
      end
      
      % Fading memory
      ws.block.previous_ssm_activeness = ws.block.reduced_ssm_activeness;
      
      reduce_ssm_data;
      
    end
    
    
    a_struct.total_step_clock(6) = a_struct.total_step_clock(6) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    % ( 3.7 ) Update the detection model
    %
    % -> done at the end only, this is why
    % we claim "online, frame-by-frame" processing
    
    if frame_ind == ws.block.stop_frame
      
      if PAR.DETECTION.DEBUG
	disp( [ '[          Updating the detection model...          ]' ] );
      end
      
      fit_model;

      a_struct.total_step_clock(7) = a_struct.total_step_clock(7) + etime( clock, tmp_clock );
      % ---------- CLOCK STEP LIMIT
      tmp_clock = clock;
      
      % Select decision thresholds

      select_thresholds;
      
    end

    a_struct.total_step_clock(8) = a_struct.total_step_clock(8) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    % ( 3.8 ) Reduce the 3-D location estimates data to a fixed size

    if PAR.DO_STC &  (  frame_ind == ws.block.stop_frame )
      
      if PAR.STC.VERBOSE
	disp( [ '[          Reducing the 3-D data to a fixed size (used by short-term clustering)...          ]' ] );
      end
      
      % Extract data samples ( delta azimuth, delta elevation )
      % for this block
      ind = find( ( sphloc_mat( 1,: ) >= ws.block.start_frame ) & ( sphloc_mat( 1,: ) <= ws.block.stop_frame ) );
      ws.block.sphloc_mat = sphloc_mat( :, ind );

      extract_stc_data;
      
      % Fading memory
      ws.block.previous_stc_data = ws.block.reduced_stc_data;
      
      % Merge with past (similar to "forgetting factor")
      % this will modify "ws.block.reduced_stc_data"
      reduce_stc_data;
      
    end
    
    a_struct.total_step_clock(9) = a_struct.total_step_clock(9) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    % ( 3.9 ) Update the short-term clustering local models
    %
    % -> done at the end only, this is why
    % we claim "online, frame-by-frame" processing

    if PAR.DO_STC  &  ( frame_ind == ws.block.stop_frame )

      if PAR.STC.VERBOSE
        disp( [ '[          Updating the short-term clustering local models...          ]' ] );
      end
    
      update_stc_model;
      
    end
    
        
    a_struct.total_step_clock(10) = a_struct.total_step_clock(10) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;


    % ----------------------------------------------------------------------
    % Sector-based MFCC extraction
    
    if PAR.DO_SBMFCC &  (  frame_ind == ws.block.stop_frame )

      % DEBUG 
      % ( 3.10 ) Sector-based binary filter in frequency domain.
      % Noise removal using USS.
      % -> preparation for next step (sector-based MFCC extraction).
      %
      % Result stored into "ws.block.abs_X_sector"
      
      if PAR.VERBOSE > 1
	disp( sprintf( [ '[          Removing noise using USS (frames %d to %d) on channel %d...          ]' ], ws.block.start_frame, ws.block.stop_frame, ws.PAR.MFCC_CHANNEL ) );
      end
      
      % Prepare for USS: reduce the data + forgetting factor
      % ( update "ws.block.datarepr" )
      update_datarepr_uss;
      
      % Apply SSM binary decision on the spectrogram,
      % for each different sector
      % ( create "ws.block.abs_X_sector" )
      extract_abs_X_sector;
      
      % Fit model on single channel spectrum,
      % and apply USS on each sector-based spectrum
      % ( modifies "ws.block.abs_X_sector" )
      remove_noise;
      
    end

    a_struct.total_step_clock( 11 ) = a_struct.total_step_clock( 11 ) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    if PAR.DO_SBMFCC &  (  frame_ind == ws.block.stop_frame )
      
      % ( 3.11 ) Extract sector-based MFCCs
      
      if PAR.VERBOSE
	disp( sprintf( [ '[          Extracting sector-based MFCCs (frames %d to %d) from filtered channel %d...     ]' ], ws.block.start_frame, ws.block.stop_frame, ws.PAR.MFCC_CHANNEL ) );
      end
      
      % Pick only the (time frame, sector) pairs that we need.
      ind = find( ( sphloc_mat( 1,: ) >= ws.block.start_frame ) & ( sphloc_mat( 1,: ) <= ws.block.stop_frame ) );

      if ~isempty( ind )
      
	nbins       = size( ws.block.abs_X_sector, 1 );
	frame_list  = sphloc_mat( 1, ind );

	block_frame_list  = repmat( frame_list - ws.block.start_frame + 1, nbins, 1 );
	sector_list       = repmat( sphloc_mat( 2, ind ), nbins, 1 );
	bin_list          = repmat( ( 1:nbins ).',        1,     numel( ind ) );
	
	tmp_ind = sub2ind( size( ws.block.abs_X_sector ), bin_list, sector_list, block_frame_list );
	% Result: tmp is a 2-D matrix ( nbins x nlocationestimates )
	tmp     = ws.block.abs_X_sector( tmp_ind );
	
	% For compatibility with "compute_mfcc",
	% give the symmetric version of the spectrum
	% The DC component is set to zero (it should not matter).
	tmp = [ zeros( 1, size( tmp, 2 )  );   tmp;  tmp( ws.PAR.FFT_FRAME_LENGTH:-1:2, : ) ];
	
	% Sanity check
	if size( tmp, 1 ) ~= 2 * ws.PAR.FFT_FRAME_LENGTH
	  error( [ mfilename ': wrong number of FFT bins.' ] );
	end
	
	% Do it!
	mfcc = compute_mfcc( tmp );

	if ~PAR.MFCC.DO_C0
	  mfcc = mfcc( 2:end, : );
	end
	
	% Sanity check
	if any( isnan( mfcc(:) ) )
	  error( [ mfilename ' is NaN-insane.' ] );
	end
	
	% Sanity check
	if any( isinf( mfcc(:) ) )
	  error( [ mfilename ' is Inf-insane.' ] );
	end
	
      	clear tmp;
	
	% Store the result
	sbmfcc_mat = [ sbmfcc_mat  [ frame_list( 1,: ); sector_list( 1,: ); mfcc ] ];

      end

    end      % if PAR.DO_SBMFCC

    % --------------------------------------------------
    % ( 3.12 ) Single-channel MFCC extraction

    if PAR.DO_MFCC  &  (  frame_ind == ws.block.stop_frame )
      
      if PAR.VERBOSE
	disp( sprintf( [ '[          Extracting single-channel MFCCs (frames %d to %d) from channel %d...            ]' ], ws.block.start_frame, ws.block.stop_frame, ws.PAR.MFCC_CHANNEL ) );
      end
      
      % Spectrum of one channel
      tmp = abs( ws.block.X( :, PAR.MFCC_CHANNEL, : ) );

      % From 3-D to 2-D
      tmp = tmp( :, : );
      
      % For compatibility with "compute_mfcc",
      % give the symmetric version of the spectrum
      % The DC component is set to zero (it should not matter).
      tmp = [ zeros( 1, size( tmp, 2 )  );   tmp;  tmp( ws.PAR.FFT_FRAME_LENGTH:-1:2, : ) ];
      
      % MFCC

      % Do it!
      mfcc = compute_mfcc( tmp( :, : ) );
      
      if ~PAR.MFCC.DO_C0
	mfcc = mfcc( 2:end, : );
      end
      
      % Sanity check
      if any( isnan( mfcc(:) ) )
	error( [ mfilename ' is NaN-insane.' ] );
      end
      
      % Sanity check
      if any( isinf( mfcc(:) ) )
	error( [ mfilename ' is Inf-insane.' ] );
      end
      
      clear tmp;
      
      % Store the result
      mfcc_mat = [ mfcc_mat  [ ws.block.start_frame:ws.block.stop_frame; mfcc ] ];
      
    end


    % ----------------------------------------------------------------------

    if PAR.VERBOSE
      chrono = chrono_check( chrono, i_frame );
    end


    a_struct.total_step_clock( 12 ) = a_struct.total_step_clock( 12 ) + etime( clock, tmp_clock );
    % ---------- STEP for CLOCK  

  end % of loop through frames

  if PAR.DO_STC
    
    % ---------- CLOCK STEP LIMIT  
    tmp_clock = clock;

    % Make sure that no sample is left behind (it may happen in the last frame)
    
    while ws.stc.n_new_sample > 0
      
      
      % Get the short-term past and future
      % This modifies: "ws.stc.past" and "ws.stc.future"
      get_stc_past_future;
      
      % Apply short-term clustering
      % This modifies the tags in "ws.stc.future"
      short_term_clustering;
      
      % Store the result
      % This modifies: "sphloc_mat", "ws.stc.last_loc_est"
      update_stc_past_future;
      
    end
    
    a_struct.total_step_clock(5) = a_struct.total_step_clock(5) + etime( clock, tmp_clock );
    % ---------- CLOCK STEP LIMIT
    
  end
  
  
  % Sanity check

  if PAR.DO_STC
    if ws.stc.n_new_sample > 0
      error( [ mfilename ': some samples were forgotten while doing the short-term clustering.' ] );
    end
  end
  
  if PAR.VERBOSE
    chrono_stop( chrono );
  end
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 4 ) Finished: return the result
  
  a_struct.sphloc_mat = sphloc_mat;
  
  if PAR.VERBOSE
    disp( [ mfilename ': done.' ] );
  end

  a_struct.total_clock     = etime( clock, start_clock );

  % Also return misc. stuff, including the last version of the local models
  
  a_struct.last_stc_model = ws.stc.model;

  a_struct.ssm = ws.ssm;
  
  if PAR.DO_SBMFCC

    % Return sector-based MFCC    
    a_struct.sbmfcc_mat = sbmfcc_mat;
    
  end

  if PAR.DO_MFCC
  
    % Return single-channel MFCC
    a_struct.mfcc_mat = mfcc_mat;
    
  end

  % And in case default parameter values were added, it is interesting to store them
  
  a_struct.PAR = ws.PAR;
  
  
% ----------------------------------------------------------------------
% Check that all mandatory parameters are available.
% Extract a few values.

function check_mandatory( a_struct, ssm, PAR )
  
  if nargin < 3
    error( [ mfilename '.check_mandatory() needs 3 input arguments.' ] );
  end
  
  required_params = { 'wavefile_list', 'start_sample', 'stop_sample' };
  check_param( required_params, fieldnames( a_struct ) );
  
  % Make it robust to slight mistakes from the user
  a_struct.start_sample = max( 1, a_struct.start_sample );
  a_size = wavread( a_struct.wavefile_list{ 1 }, 'size' );
  a_struct.stop_sample  = min( a_size( 1 ), a_struct.stop_sample );
  
  required_params = { 'framelength', 'fs', 'c', 'pairs', 'Z', 'w', 'geometry', 'sector' };
  check_param( required_params, fieldnames( ssm ) );
  
  check_param( { 'init_sph' }, fieldnames( ssm.sector ) );
  
  required_params = { 'FRAMELENGTH_SEC', 'TIMEFRAME_OVERLAP', 'DETECTION', 'BLOCK_SIZE_SEC', 'SPECTRUM_DO_ZERO_MEAN', 'SPECTRUM_PREEMP', 'FASTTDE_SQUARE_SUBARRAY_LIST', 'FASTTDE_UPFACTOR', 'DO_SBMFCC', 'DO_MFCC'  };
  check_param( required_params, fieldnames( PAR ) );
    
  % Check subarrays' definitions
  if ~iscell( PAR.FASTTDE_SQUARE_SUBARRAY_LIST )
    error( [ mfilename '.check_mandatory(): PAR.FASTTDE_SQUARE_SUBARRAY_LIST must be a cell array.' ] );
  end
  
  nchannel  = numel( a_struct.wavefile_list );
  
  nsubarray = numel( PAR.FASTTDE_SQUARE_SUBARRAY_LIST );
  
  for j = 1:nsubarray
  
    siz = size( PAR.FASTTDE_SQUARE_SUBARRAY_LIST{ j } );
    do_error = 0;
    if numel( siz ) ~= 2
      do_error = 1;
    else
      if ~all( siz == [ 2 2 ] );
	do_error = 1;
      end
    end
    if do_error
      error( [ mfilename '.check_mandatory(): PAR.FASTTDE_SQUARE_SUBARRAY_LIST must be a cell array of 2 x 2 matrices.' ] );
    end
    
    if ~all( ismember( PAR.FASTTDE_SQUARE_SUBARRAY_LIST{ j }(:), 1:nchannel ) )
      error( [ mfilename '.check_mandatory(): PAR.FASTTDE_SQUARE_SUBARRAY_LIST contents must all be between 1 and nchannel.' ] );
    end
    
  end


  
  if PAR.DO_MFCC | PAR.DO_SBMFCC
    
    check_param( { 'MFCC_CHANNEL' }, fieldnames( PAR ) );
    
    if ( PAR.MFCC_CHANNEL < 1 ) | ( PAR.MFCC_CHANNEL > numel( a_struct.wavefile_list ) )
      error( [ mfilename '.check_mandatory(): "PAR.MFCC_CHANNEL" must be in [ 1 nchannels ].' ] );
    end
    
  end
  
  check_param( { 'FAR_T', 'FRR_T', 'SPACE_TIME_DILATION' }, fieldnames( PAR.DETECTION ) );
  
  required_params = { 'STC' };
  check_param( required_params, fieldnames( PAR ) );
  
  required_params = { 'T_SHORT' };
  check_param( required_params, fieldnames( PAR.STC ) );
  
  % Basic consistency check
  
  for i_channel = 1:numel( a_struct.wavefile_list )
    
    [ tmp, fs ] = wavread( a_struct.wavefile_list{ i_channel }, 'size' );
    
    if fs ~= ssm.fs
      error( [ mfilename '.check_mandatory(): inconsistent sampling frequencies.' ] );
    end
    
  end
  
  
% ----------------------------------------------------------------------
% Check for optional parameters.
% Fill in default values.

function PAR = check_optional( PAR, ssm )
  
  if nargin < 2
    error( [ mfilename '.check_optional() needs 1 input argument.' ] );
  end
  
  
  PAR_DEFAULT.FFT_FRAME_LENGTH = ssm.framelength;

  % Derive framelength and frameshift in seconds
    
  PAR.FRAMESHIFT_SEC  = PAR.FRAMELENGTH_SEC * ( 1 - PAR.TIMEFRAME_OVERLAP );
  PAR.FRAMESHIFT_SEC  = round( PAR.FRAMESHIFT_SEC * ssm.fs ) / ssm.fs;
  
  
  % ( 1 ) Sector-based detection

  DETECTION_DEFAULT.MIN_ITER = 5;
  DETECTION_DEFAULT.MAX_ITER = 20;

  DETECTION_DETAULT.N_SECTORS_MAX = +Inf;
  
  % Data representation: give equal weight to:
  % - all data before the current block.
  % - all data in the current block.
  DETECTION_DEFAULT.DATAREPR_MAXPREVIOUS_SEC       = PAR.BLOCK_SIZE_SEC;  % In seconds
  DETECTION_DEFAULT.DATAREPR_MAXCURRENT_SEC        = PAR.BLOCK_SIZE_SEC;  % In seconds
  
  % Parameters for the EM fitting
  opt = optimset( @fminbnd );
  opt = optimset( opt, 'Display', 'none', 'MaxFunEvals', 60, 'MaxIter', 30, 'TolX', +Inf, 'TolFun', 1e-4 );
  DETECTION_DEFAULT.FMINSEARCH_OPTIONS = opt;
  
  DETECTION_DEFAULT.AUTO_INIT_N = 30;
  DETECTION_DEFAULT.LKLD_CONV   = 1e-10;
  DETECTION_DEFAULT.WD_WG_FLOOR = 1e-10;
  DETECTION_DEFAULT.A0_A1_FLOOR = 1e-10;
  DETECTION_DEFAULT.B0_B1_FLOOR = 1e-10;
  
  
  
  % Verbosity, mostly to see how the EM converges
  DETECTION_DEFAULT.VERBOSE = 1;
  
  % Mostly to show on a figure the convergence of the models after EM
  DETECTION_DEFAULT.DEBUG   = 0;
  DETECTION_DEFAULT.DEBUG_HIST_NBINS = 100;
  
  if ~isfield( PAR, 'DETECTION' )
    PAR.DETECTION = DETECTION_DEFAULT;
  else
    PAR.DETECTION = fill_default( PAR.DETECTION, DETECTION_DEFAULT );
  end
  

  % ( 2 ) GCC-PHAT Time-Delay Estimation

  % To define the length of the FIR filter used for upsampling the GCC-PHAT
  
  PAR_DEFAULT.FASTTDE_FILTERORDER = 6 * PAR.FASTTDE_UPFACTOR;
  

  % ( 4 ) Unsupervised Spectral Subtraction (USS)
  
  USS_DEFAULT = [];

  % Type of model: the Shifted Erlang starts at
  %   ERLANG_T_FACTOR * sigma_{rayleigh}
  %
  % Must match the name of this script
  USS_DEFAULT.ERLANG_T_FACTOR = 1.0;
  
  % For magnitude data representation (for USS):
  % we'll take some well-chosen data points from the past (previous blocks),
  % some well-chosen data points from the future (next block),
  % and fit the Rayleigh + Shifted Erlang model on the whole.
  %
  % The small number of points (100 typically) allows for quick EM fitting.
  % It does not hurt since the model is simple (bi-modal in 1-D).
  %
  % >> Rule of thumb is: keep both values equal. <<
  
  USS_DEFAULT.data_par.past_npoints_max = 100;
  USS_DEFAULT.data_par.future_npoints   = USS_DEFAULT.data_par.past_npoints_max;
  
  % Parameters for EM fitting
  
  USS_DEFAULT.em_par.erlang_t_factor = USS_DEFAULT.ERLANG_T_FACTOR;
  USS_DEFAULT.em_par.min_iter        = 5;
  USS_DEFAULT.em_par.max_iter        = 200;
  USS_DEFAULT.em_par.sigma_floor     = 1e-10;
  USS_DEFAULT.em_par.lambda_floor    = 1e-10;
  
  %   "moment" if do_fminsearch is 0  -> faster and better in noise
  %   "ml" if do_fminsearch is 1      -> exact but worse in noise
  
  USS_DEFAULT.em_par.do_fminsearch = 0;
  
  opt = optimset( @fminsearch );
  opt = optimset( opt, 'Display', 'none', 'MaxFunEvals', 50, 'MaxIter', 25, 'TolX', 1e-4, 'TolFun', 1e-4 );
  USS_DEFAULT.em_par.fminsearch_options = opt;    

  % Verbosity (especially iterations during training
  USS_DEFAULT.em_par.verbose = 0;
  
  % Check for NaNs of Infs (this particular one does not cost much)
  USS_DEFAULT.em_par.sanity_check = 1;

  if ~isfield( PAR, 'USS' )
    PAR.USS = USS_DEFAULT;
  else
    PAR.USS = fill_default( PAR.USS, USS_DEFAULT );
  end
  
  % ( 5 ) MFCCs
  
  MFCC_DEFAULT = [];

  % Check for NaNs or Infs, set this to zero when you are confident
  MFCC_DEFAULT.sanity_check = 0;
  
  MFCC_DEFAULT.B = 24;   % Number of MEL filter banks

  % Number of non-zero cepstrum coefficients (C1,C2,...,CN)
  % Default = maximum = B-1
  MFCC_DEFAULT.N = MFCC_DEFAULT.B-1;   

  % Cepstral liftering: 0 means none, [] means default value
  MFCC_DEFAULT.L = [];   
  
  MFCC_DEFAULT.DO_C0 = 1;  % Do we want C0?

  % ( normally leave both to zero )
  MFCC_DEFAULT.usepower = 0;  % Energy (0) or power (1)
  MFCC_DEFAULT.verbose  = 0;

  % Flooring is not needed given that we do USS (see "function remove_noise" below)
  MFCC_DEFAULT.fbank_floor = 0;

  if ~isfield( PAR, 'MFCC' )
    PAR.MFCC = MFCC_DEFAULT;
  else
    PAR.MFCC = fill_default( PAR.MFCC, MFCC_DEFAULT );
  end

  % Default value for liftering
  if isempty( PAR.MFCC.L )
    PAR.MFCC.L = 2 * ( PAR.MFCC.N - 1 );
  end
  
  
  % ( 6 ) Verbosity at my level
  
  PAR_DEFAULT.VERBOSE = 1;

  % ( 6 ) Short-term clustering
  
  PAR_DEFAULT.DO_STC  = 1;  % by default, do it!
  
  % Fill in where missing
  PAR = fill_default( PAR, PAR_DEFAULT );
  
  % FASTTDE: check that the filter order is even
  if mod( PAR.FASTTDE_FILTERORDER, 2 )
    error( [ mfilename '.check_optional(): "PAR.FASTTDE_FILTERORDER" must be even.' ] );
  end
  
  % Default parameters for short-term clustering
  if PAR.DO_STC
    
    STC_DEFAULT.VERBOSE = 0;
    
    % Order of magnitude of the maximum acceptable angle difference, in radians
    % -> this parameters is used only to process the first block of data.
    STC_DEFAULT.MAX_DAZ =  5/180*pi;
    
    
    STC_DEFAULT.N_PAST_FRAMES    = [];  % default value is "PAR.STC.T_SHORT"

    STC_DEFAULT.N_FUTURE_SAMPLES = 6;
    
    % To prune impossible merges between past clusters and future clusters:
    % any graph with a local same/not sure  decision less likely than this
    % threshold is pruned out.
    %    
    % -Inf is not advised since the number of samples in the past may be large,
    % unless you are sure that you have low complexity (e.g. N_FUTURE_SAMPLES = 1).
    STC_DEFAULT.LLR_THRESHOLD = -20;

    % Graphs: full search by default
    STC_DEFAULT.N_GRAPHS_MAX  = +Inf;
    STC_DEFAULT.N_MERGES_MAX  = +Inf;
    
    STC_DEFAULT.MODEL_STD_FLOOR = 1e-3;  % Floor on the Gaussian std dev, in radians
    
    % For initialization when updating the local models
    % For EM training: set some small value to prevent "cyclic" behaviours
    % being captured as "same" (e.g. PAR.STC.DAZ_MAX = 5 / 180 * pi for a maximum of 5 degrees).
    STC_DEFAULT.DAZ_MAX = +Inf;
    
    STC_DEFAULT.AUTO_INIT_N   = 10;
    
    PAR.STC = fill_default( PAR.STC, STC_DEFAULT );
    
    
    if isempty( PAR.STC.N_PAST_FRAMES )
      PAR.STC.N_PAST_FRAMES   = PAR.STC.T_SHORT;
    end
    
    
    if ( PAR.STC.N_GRAPHS_MAX == +Inf ) & ( PAR.STC.N_FUTURE_SAMPLES > 11 )
      disp( [ mfilename ' short-term clustering will never end. You may want to limit the number of graphs.' ] );
    end
    
  end

% ----------------------------------------------------------------------
% Get details such as framelength etc.

function  a_struct = get_details( a_struct, ssm, PAR )

  if nargin < 3
    error( [ mfilename '.get_details() needs 3 input arguments.' ] );
  end

  a_struct.framelength       = PAR.FRAMELENGTH_SEC * ssm.fs;  % In samples
  a_struct.framelength_sec   = PAR.FRAMELENGTH_SEC;   % In seconds
  
  % FFT might be different (e.g. to reduce computational cost)
  a_struct.fft_framelength   = ssm.framelength;
  
  a_struct.frameshift_sec    = a_struct.framelength_sec * ( 1 - PAR.TIMEFRAME_OVERLAP );  % In seconds
  a_struct.frameshift        = round( a_struct.frameshift_sec * ssm.fs );  % In samples
  a_struct.frameshift_sec    = a_struct.frameshift / ssm.fs;

  a_struct.block_size_frames      = round( PAR.BLOCK_SIZE_SEC       / a_struct.frameshift_sec );
  a_struct.init_block_size_frames = round( PAR.INIT_BLOCK_SIZE_SEC  / a_struct.frameshift_sec );
  
  a_struct.start_frame       = 1 + ceil( ( a_struct.start_sample - 1 ) / a_struct.frameshift );
  a_struct.stop_frame        = 1 + floor( ( a_struct.stop_sample - a_struct.framelength )  / a_struct.frameshift );
  a_struct.n_frames          = 1 + a_struct.stop_frame - a_struct.start_frame;
  
  a_struct.n_channels        = numel( a_struct.wavefile_list );
  a_struct.win               = diag( sparse( hamming( a_struct.framelength, 'periodic' ) ) );
  a_struct.n_pairs           = size( ssm.pairs, 1 );
  a_struct.n_sectors         = numel( ssm.sector );
  
  a_struct.conj_ssm_Z        = conj( ssm.Z );


  % Define how we store past data, for detection
  
  a_struct.datarepr_maxprevious_frames = max( 1, round( PAR.DETECTION.DATAREPR_MAXPREVIOUS_SEC / a_struct.frameshift_sec ) );
  
  a_struct.datarepr_maxcurrent_frames  = max( 1, round( PAR.DETECTION.DATAREPR_MAXCURRENT_SEC / a_struct.frameshift_sec ) );

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % FASTTDE: convert each subarray into a list of two pairs,
  % and create the list of pairs

  a_struct.fasttde.square_subarray_i_pair_list = cell( size( PAR.FASTTDE_SQUARE_SUBARRAY_LIST ) );

  a_struct.fasttde.nsubarray = numel( PAR.FASTTDE_SQUARE_SUBARRAY_LIST );
  
  for j = 1:a_struct.fasttde.nsubarray
    
    % sort columns
    tmp = sort( PAR.FASTTDE_SQUARE_SUBARRAY_LIST{ j }.' ).';
    
    % Find the index of the first pair, in "ssm.pairs"

    i_pair_1 = find( ( ssm.pairs(:,1) == tmp(1,1) ) & ( ssm.pairs(:,2) == tmp(1,2) ) );
    
    if numel( i_pair_1 ) ~= 1
      error( [ mfilename '.get_details(): sorry! Could not find this pair within the list predefined for SAM-SPARSE-MEAN.' ] );
    end
    
    % Find the index of the second pair, in "ssm.pairs"

    i_pair_2 = find( ( ssm.pairs(:,1) == tmp(2,1) ) & ( ssm.pairs(:,2) == tmp(2,2) ) );
    
    if numel( i_pair_2 ) ~= 1
      error( [ mfilename '.get_details(): sorry! Could not find this pair within the list predefined for SAM-SPARSE-MEAN.' ] );
    end
    
    % Store the two indices
    % This will be very useful to get the part of the GCCPHAT that we are interested in

    a_struct.fasttde.square_subarray_i_pair_list{ j } = [ i_pair_1  i_pair_2 ];
    
  end
  
  % Stack all the pairs we need
  % This will be also useful to do inverse FFT on what we need only
  
  a_struct.fasttde.i_pair_list = [ a_struct.fasttde.square_subarray_i_pair_list{:} ];
  a_struct.fasttde.n_pair      = numel( a_struct.fasttde.i_pair_list );

  % Finally, determine where in the inverse FFT we'll look at,
  % for each subarray. Each value is an index in "a_struct.fasttde.i_pair_list".
  
  a_struct.fasttde.i_i_pair_list = zeros( a_struct.fasttde.nsubarray, 2 );
  
  for j = 1:a_struct.fasttde.nsubarray
    
    i_pair = a_struct.fasttde.square_subarray_i_pair_list{ j };
    
    a_struct.fasttde.i_i_pair_list( j, : ) = [ find( i_pair( 1 ) == a_struct.fasttde.i_pair_list ), find( i_pair( 2 ) == a_struct.fasttde.i_pair_list ) ];
    
  end
  
  
  % For each subarray, compute the homogenous (4-row) matrix of rigid transformation (translation + rotation)
  
  a_struct.fasttde.square_subarray_Pmat_list = cell( size( PAR.FASTTDE_SQUARE_SUBARRAY_LIST ) );
  
  a_struct.fasttde.square_subarray_list      = PAR.FASTTDE_SQUARE_SUBARRAY_LIST;
  
  for j = 1:a_struct.fasttde.nsubarray
    
    sqsub = a_struct.fasttde.square_subarray_list{ j };
    
    % Unit vectors defining the local coordinate system, as in Brandstein's thesis
    
    x_j = diff( ssm.geometry( :, sqsub( 1, [2 1] ) ), [], 2 );
    x_j = x_j / sqrt( dot( x_j, x_j ) );
    
    y_j = diff( ssm.geometry( :, sqsub( 2, [2 1] ) ), [], 2 );
    y_j = y_j / sqrt( dot( y_j, y_j ) );
    
    % Sanity check
    if abs( dot( x_j, y_j ) > 100 * eps )
      error( [ mfilename ': each subarray must be composed of two ORTHOGONAL microphone pairs.' ] );
    end
    
    z_j = cross( x_j, y_j );
    z_j = z_j / sqrt( dot( z_j, z_j ) );
    
    % Translation w.r.t. the original coordinate systen
    
    m_j = mean( ssm.geometry( :, sqsub(:).' ), 2 );
    
    % Store both rotation and translation in a 4-row matrix of homogenous coordinates
    
    Pmat = [ [ x_j(:)  y_j(:)  z_j(:)  m_j(:) ]; 0 0 0 1 ];
    Pmat( find( abs( Pmat ) < 100 * eps ) ) = 0;
    
    a_struct.fasttde.square_subarray_Pmat_list{ j } = sparse( Pmat );
    
  end


  % FASTTDE: Prepare low-pass filter for upsampling.
  % This filter is an approximated sinc filter.
  % Note that we will only need to interpolate
  % a small portion of the time-domain GCC-PHAT,
  % centered around time-delay zero and bounded
  % by a maximum time-delay corresponding to 
  % the distance between the two microphones of each pair.

  a_struct.fasttde.upfactor = PAR.FASTTDE_UPFACTOR;
  a_struct.fasttde.filterorder = PAR.FASTTDE_FILTERORDER;

  a_struct.fasttde.B         = PAR.FASTTDE_UPFACTOR * fir1( PAR.FASTTDE_FILTERORDER,  1 / PAR.FASTTDE_UPFACTOR );
  a_struct.fasttde.halforder = floor( PAR.FASTTDE_FILTERORDER / 2 );
  
  % Determine the portion of time-domain GCC-PHAT
  % that we will look at
  
  a_struct.fasttde.maxdelay = ssm.fs / ssm.c * pair_distance( ssm.geometry, ssm.pairs( a_struct.fasttde.i_pair_list, : ) );
  
  a_struct.fasttde.halfportion = 1 + ceil( max( a_struct.fasttde.maxdelay ) ) + ceil( PAR.FASTTDE_FILTERORDER / PAR.FASTTDE_UPFACTOR );
 
  if a_struct.fasttde.halfportion >= a_struct.framelength
    error( [ mfilename '.get_details(): maximum time-delay too large, please check "ssm.geometry" or increase "PAR.FRAMELENGTH_SEC".' ] );
  end
  
  % In the original time-domain GCC-PHAT
  % ( in the 2D matrix called "gccphat" further below )
  
  a_struct.fasttde.rowzero  = 1 + a_struct.framelength;   % Where time-delay zero is in the "gccphat" matrix
  a_struct.fasttde.rowstart = a_struct.fasttde.rowzero - a_struct.fasttde.halfportion;
  a_struct.fasttde.rowend   = a_struct.fasttde.rowzero + a_struct.fasttde.halfportion;
  a_struct.fasttde.nrows    = 1 + a_struct.fasttde.rowend - a_struct.fasttde.rowstart;
  
  % After upsampling
  % ( in the 2D matrix called "up_gccphat" further below )
  
  a_struct.fasttde.up_rowzero = a_struct.fasttde.halforder + 1 + ( a_struct.fasttde.rowzero - a_struct.fasttde.rowstart ) * PAR.FASTTDE_UPFACTOR;
  a_struct.fasttde.up_nrows   = a_struct.fasttde.nrows * PAR.FASTTDE_UPFACTOR;
  
  % Corresponding time-delays
  % -> same convention as in "location_to_timedelay.m":
  %    within each pair,
  %    these are "advance" of mic. #1 w.r.t. mic. #2
  a_struct.fasttde.up_row_td = ( a_struct.fasttde.up_rowzero - (1:a_struct.fasttde.up_nrows) ) / PAR.FASTTDE_UPFACTOR;
  

  % For each sector and each microphone pair,
  % store minimum and maximum time-delays in samples

  nsectors = numel( ssm.sector );
  a_struct.fasttde.td_min = repmat( NaN, a_struct.fasttde.n_pair, nsectors );
  a_struct.fasttde.td_max = repmat( NaN, a_struct.fasttde.n_pair, nsectors );

  a_struct.fasttde.up_range_def = repmat( -Inf, [ a_struct.fasttde.up_nrows, a_struct.fasttde.n_pair, nsectors ] );
  
  
  for a = 1:nsectors
  
    s = ssm.sector( a );

    % Create a grid of locations covering this sector
    az_list = s.az_min + ( s.az_max - s.az_min ) * (0:20) / 20;
    el_list = s.el_min + ( s.el_max - s.el_min ) * (0:20) / 20;
    r_list  = exp( log( s.r_min ) + ( log( s.r_max / s.r_min ) * (0:20) / 20 ) );

    % Move to cartesian
    x = sph2cart_mat( peigne_vec( az_list, el_list, r_list ) );
  
    % Move to time-delays
    td = location_to_timedelay( ssm.geometry, ssm.pairs( a_struct.fasttde.i_pair_list, : ), x, ssm.fs, ssm.c );
  
    a_struct.fasttde.td_min( :, a ) = min( td, [], 2 );
    a_struct.fasttde.td_max( :, a ) = max( td, [], 2 );
    
  end
  
  % For each sector and each microphone pair,
  % convert the range of time-delays [ td_min, td_max ] 
  % from sample values to row indices in the upsampled time-domain GCCPHAT matrix 
  
  a = round( a_struct.fasttde.up_rowzero - a_struct.fasttde.td_min * a_struct.fasttde.upfactor );
  b = round( a_struct.fasttde.up_rowzero - a_struct.fasttde.td_max * a_struct.fasttde.upfactor );
  
  a_struct.fasttde.up_index_min = min( a, b );
  a_struct.fasttde.up_index_max = max( a, b );
  
  % Convert each pair of value [ up_index_min( i,j )   up_index_max( i,j ) ] into a vector
  % of zeroes and "-Inf"s that will be used to "mask" the time-domain GCCPHAT matrix
  % in the "one_array_fasttde_dirloc" function.
  
  a_struct.fasttde.up_range_def = repmat( -Inf, [ a_struct.fasttde.up_nrows, a_struct.fasttde.n_pair, nsectors ] );
  
  for i_sec = 1:nsectors
    
    for i_pair = 1:a_struct.fasttde.n_pair
      
      ind = a_struct.fasttde.up_index_min( i_pair, i_sec ):a_struct.fasttde.up_index_max( i_pair, i_sec );
      a_struct.fasttde.up_range_def( ind, i_pair, i_sec ) = 0;
      
    end
    
  end
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Parameters for MFCC extraction 
  % (both sector-based MFCCs and single-channel MFCCs)
  
  if PAR.DO_MFCC | PAR.DO_SBMFCC

    % calculate coefficients for the filter bank analysis
    % ( magnitude spectrogram -> MEL fbank energy spectrogram )
    
    a_par = [];
    a_par.fs               = ssm.fs;
    a_par.fft_framelength  = PAR.FFT_FRAME_LENGTH;
    a_par.mfcc_band_hz     = [ 50  +Inf ];   % 50 Hz high-pass
    a_par.mfcc_B           = PAR.MFCC.B;
    
    a_struct.fbank_mat     = prepare_fbank_mat( a_par );
    
    % calculate coefficients for the DCT analysis
    % ( log MEL fbank energy spectrogram -> MFCC coefficients )

    a_par = [];
    a_par.mfcc_B       = PAR.MFCC.B;
    a_par.mfcc_N       = PAR.MFCC.N;
    a_par.mfcc_L       = PAR.MFCC.L;

    a_struct.dct_mat   = prepare_dct_mat( a_par );
    
  end

  
% ----------------------------------------------------------------------
% Find all possible connected graphs of N points,
% where each possible edge belongs to the  "same" or "different" category.
% This is a recursive implementation: it calls itself with N-1.
% The difficulty is to make sure that each graph is consistent
% ( if A = B and B = C then A = C ).
%
% Output: GRAPH_LIST_N contains a list of graphs with N points,
% stored as a cell array of N - by - N matrices of binary values
% In each matrix, only the upper triangular part is used,
% where a 1 indicates "same" and a 0 indicates "different".
%
% An optional pruning is possible, through matrices "permitted_0" and
% "permitted_1". Each is a N - by - N matrix of ones (possible link)
% and zeros (forbidden link). Hence zeroes form the prior knowledge
% (default: ones(N,N)). Only the upper triangular part is used.
%
% N_GRAPHS_MAX is the maximum number of graphs to construct.
% If reached, we issue a warning and return the N_GRAPHS_MAX graphs.
% default value: +Inf.
%
% Each combination is stored as a sparse matrix ( only strictly upper triangular
% part is used, lower triangular part is left as zeroes ).

function graph_list_n = construct_graph_list_pruned( n, graph_list_n_minus_1, permitted_0, permitted_1, n_graphs_max )
  
  if nargin < 1
    error( sprintf( [ mfilename '.construct_graph_list_pruned[%d] needs at least one parameter' ] ,n) );
  end
  
  if ~isnumeric( n )
    error( sprintf( [ mfilename '.construct_graph_list_pruned[%d] needs an integer N >= 0' ] ,n) );
  end
  
  if rem( n, 1 ) | ( n < 0 )
    error( sprintf( [ mfilename '.construct_graph_list_pruned[%d] needs an integer N >= 0' ] ,n) );
  end
  
  n = floor( n );
  
  % No point: empty matrix
  if n < 1
    graph_list_n = { sparse( [], [], [], 0, 0 ) };
    return;
  end
  
  % List of graphs on which we build on
  % Default: empty
  if ~exist( 'graph_list_n_minus_1', 'var' )
    graph_list_n_minus_1 = [];
  end
  
  % Limitation of the permitted "different" edges
  % permitted_0 is a N x N matrix of binary values,
  % where "1" on an edge indicates a permitted "different" decision for this edge,
  % and   "0" on an edge indicates a forbidden "different" decision for this edge,
  if ~exist( 'permitted_0', 'var' )
    permitted_0 = [];
  end
  
  % Limitation of the permitted "same" edges
  % permitted_1 is a N x N matrix of binary values,
  % where "1" on an edge indicates a permitted "same" decision for this edge,
  % and   "0" on an edge indicates a forbidden "same" decision for this edge,
  if ~exist( 'permitted_1', 'var' )
    permitted_1 = [];
  end
  
  if ~exist( 'n_graphs_max', 'var' )
    n_graphs_max = [];
  end
    
  %%
  
  if isempty( permitted_0 )
    permitted_0 = ones( n );
  end
  
  if isempty( permitted_1 )
    permitted_1 = ones( n );
  end
  
  if isempty( n_graphs_max )
    n_graphs_max = +Inf;
  end
  
  % 1 or more points:
  % first make sure we have the result from the previous order
  if isempty( graph_list_n_minus_1 )
    graph_list_n_minus_1 = construct_graph_list_pruned( n - 1, [], permitted_0(1:n-1, 1:n-1), permitted_1(1:n-1, 1:n-1), n_graphs_max );
  end
  
  % Quick sanity check
  if size( graph_list_n_minus_1{ 1 }, 1 ) ~= n-1
    error( [ mfilename '.construct_graph_list_pruned() is insane' ] );
  end
  
  % Make sure we have not reached the maximum number of graphs yet
  if length( graph_list_n_minus_1 ) >= n_graphs_max
    disp( sprintf( 'construct_graph_list_pruned[%d]: warning! graph_list_n_minus_1 already has reached maximum number of graphs.', n ) );
    graph_list_n = graph_list_n_minus_1;
    return;
  end
  
  % Algorithm to construct all possible graphs of N points:
  %
  % For each possible graph of N-1 points:
  %
  %   Find the G PERMITTED cliques of connected points.
  %   PERMITTED cliques means that all its vertices of the clique have ones 
  %   in PERMITTED_0_POINTS and in PERMITTED_1_POINTS.
  %   There are G possible ways to add the nth point to an existing clique.
  
  permitted_0_points = permitted_0( 1:n-1, n );
  permitted_1_points = permitted_1( 1:n-1, n );
  
  graph_list_n = {};
  
  for a = 1:length( graph_list_n_minus_1 )
    
    prev_combi  = graph_list_n_minus_1{ a };
    [ link_i, link_j ] = find( prev_combi );
    
    point_visited = zeros( 1, n - 1 );
    
    cliques = {};
    
    
    % 1) Build cliques of connected points
    % for this combination
    current_clique = [];
    while ~all( point_visited )

      current_clique_size = length( current_clique );
      
      if isempty( current_clique )
	
	% Create new clique
	aux = find( ~point_visited );
	current_clique = [ aux( 1 ) ];
	point_visited( aux( 1 ) ) = 1;
	
      else
	
	% Expand current clique
	v = [];
	v = [ v link_j( find( ismember( link_i, current_clique ) ) ).' ];
	v = [ v link_i( find( ismember( link_j, current_clique ) ) ).' ];
	if ~isempty( v )
	  current_clique = unique( [ current_clique, v ] );
	  point_visited( v ) = 1;
	end
	
      end
      
      % Current clique finished?
      if ( length( current_clique ) == current_clique_size ) | all( point_visited )
	% Yes, store it and start a new one
	cliques{ end+1 } = current_clique;
	current_clique = [];
      end
      
    end
    
    % 2) Which cliques are permitted?
    
    permitted_cliques = [];
    for b = 1:length( cliques )

      if all( permitted_1( cliques{ b } ) ) & all( permitted_0( ~cliques{ b } ) )
	permitted_cliques = [ permitted_cliques  b ];
      end
    
    end
    
    % 3) Based on these permitted cliques build new combinations
    
    G = length( permitted_cliques );

    % Set the right size
    prev_combi( n, n ) = 0;
    
    % Don't forget the case where the point is isolated
    graph_list_n{ end + 1 } = prev_combi;
    
    for b = 1:G
      
      the_clique = cliques{ permitted_cliques( b ) };
      p = length( the_clique );
      graph_list_n{ end + 1 } = prev_combi + sparse( the_clique, n * ones( 1, p ), ones( 1, p ), n, n );
    
      if length( graph_list_n ) > n_graphs_max
	disp( sprintf( ['construct_graph_list_pruned[%d]: warning! reached maximum number of graphs.'], n) );
	graph_list_n = { graph_list_n{ 1:n_graphs_max } };
	return;
      end
      
    end
    
  end
  

% ----------------------------------------------------------------------
% Get the spectrum, then compute the complex-domain GCC-PHAT

function get_observed_gccphat
  
  global ws
  
  start_frame = ws.block.start_frame;
  stop_frame  = ws.block.stop_frame;
  
  n_frames = 1 + stop_frame - start_frame;
  
  ws.block.n_frames = n_frames;
  
  % Load samples
  
  start_sample = ( start_frame - 1 ) * ws.a_struct.frameshift;  % may be zero ( additional sample for pre-emphasis )
  stop_sample  = ws.a_struct.framelength + ( stop_frame - 1 ) * ws.a_struct.frameshift;
  n_samples    = stop_sample - start_sample + 1;
  
  % Where we'll store the FFT of each channel
  % We'll store strictly positive frequencies only
  
  ws.block.X = zeros( [ ws.a_struct.fft_framelength, ws.a_struct.n_channels, n_frames ] );
  
  for i_channel = 1:ws.a_struct.n_channels
    
    sig = wavread( ws.a_struct.wavefile_list{ i_channel }, [ max( 1, start_sample )  stop_sample ] );
    
    if numel( sig ) < n_samples
      sig = [ 0; sig ];
    end
    
    % Bufferize
    
    sig = buffer( sig, 1 + ws.a_struct.framelength, 1 + ws.a_struct.framelength - ws.a_struct.frameshift, 'nodelay' );
    
    % Zero-mean
    
    if ws.PAR.SPECTRUM_DO_ZERO_MEAN
      sig = sig - repmat( mean( sig, 1 ), size( sig, 1 ), 1 );
    end
    
    % Pre-emphasis
    sig = sig( 2:end, : ) - ws.PAR.SPECTRUM_PREEMP * sig( 1:end-1, : );
    
    % zero-padded FFT of hamming-windowed signal
    tmp = fft( ws.a_struct.win * sig, 2 * ws.a_struct.fft_framelength, 1 );
    ws.block.X( :,i_channel,: ) = permute( tmp( 2:ws.a_struct.fft_framelength+1, : ), [ 1 3 2 ] );
    
  end
  
  clear sig tmp

  
  % From FFT to GCC-PHAT
  % Knapp-Carter 1976
  ws.block.ssm_observed_gccphat = complex( zeros( [ ws.a_struct.fft_framelength  ws.a_struct.n_pairs  n_frames ] ) );


  if 1

    % C implementation (about 26 times faster, and less memory-consuming)
    compute_gccphat_c;

  else

    % Original MATLAB implementation
    % Kept here to summarize what the C implementation does
    ws.block.ssm_observed_gccphat = exp( j * angle( ws.block.X( :, ws.ssm.pairs( :,1 ), : ) .* conj( ws.block.X( :, ws.ssm.pairs( :,2 ), : ) ) ) );

  end
  
  % Keep "X" in memory, it will be used later
  
  
% ----------------------------------------------------------------------
% Get SSM activeness from complex domain GCC-PHAT
% measurements, store it into "ws.block.ssm_activeness"
% ( as well as the underlying binary decision
% in "ws.block.ssm_binary_decision" ).
%
% SSM stands for SAM-SPARSE-MEAN (IDIAP RR 04-67).

function get_ssm_activeness
  
  global ws
  
  n_frames = size( ws.block.ssm_observed_gccphat, 3 );
  
  
  ws.block.ssm_binary_decision = uint8( zeros( ws.a_struct.fft_framelength, ws.a_struct.n_sectors, n_frames ) );
  ws.block.ssm_activeness      = zeros( ws.a_struct.n_sectors, n_frames );
  
  if 1

    % C implementation (20 times faster)
    %
    % Note that the outputs must be allocated, as was done a few lines above.

    get_ssm_activeness_c;

  else


    % Original MATLAB implementation
    % Kept here to summarize what the C implementation does.

    % To do it in one step, we would need a large amount of memory
    % ( 4-D )
    % -> loop through frames
    
    for i_frame = 1:n_frames
      
      % "MEAN": Average delay-sum power in each sector,
      % at each frequency
      
      adsp = permute( mean( 0.5 * ( 1 + real( repmat( ws.block.ssm_observed_gccphat( :,:,i_frame ), [ 1 1 ws.a_struct.n_sectors ] ) .* ws.a_struct.conj_ssm_Z ) ), 2 ), [ 1 3 2 ] );
      
      % "SPARSE": Sparsity assumption,
      % in each frequency bin, only the dominant sector is active.
      
      binary_decision = ( adsp == repmat( max( adsp, [], 2 ), 1, ws.a_struct.n_sectors ) );
      
      % If two sectors have same "adsp" value for a frequency bin, 
      % we ignore that frequency ( rare case ).
      % This way we still follow the sparsity constraint
      % ( at most one sector active in a given bin ).
      
      ind_violation = find( sum( binary_decision, 2 ) > 1 );
      binary_decision( ind_violation, : ) = 0;
      
      % -- store the binary decision
      ws.block.ssm_binary_decision( :,:, i_frame ) = binary_decision;
      
      % -- store the activeness = bandwidth occupied by each sector
      ws.block.ssm_activeness( :, i_frame ) =  mean( binary_decision, 1 ).';
      
    end

  end

  
% ----------------------------------------------------------------------
% Take the detection decision: 0 or 1, for each (sector, time frame)
%
% Input:
% - ws.block.ssm_log_post_act
% - ws.previous_block_ssm_log_post_act ( if it exists )
% - ws.block.far_threshold
% - ws.block.frr_threshold
% - ws.PAR.DETECTION.SPACE_TIME_DILATION
%
% Output:
% - ws.ssm_detection_decision( :, ws.block.start_frame:ws.block.stop_frame )

function take_detection_decision
  
  global ws
  
  % Get the last few frames, if any
  
  tmp_lpa = [];
  if isfield( ws, 'previous_block_ssm_log_post_act' )
    tmp_lpa = ws.previous_block_ssm_log_post_act;
  else
    ws.previous_block_ssm_log_post_act = [];
  end
  
  % Get the frames of the current block
  tmp_lpa = [ tmp_lpa  ws.block.ssm_log_post_act ];
  
  
  
  % ( 1 ) Consider the most conservative threshold  (FAR = FAR_T)
  % Idea: in each utterance of speech, we assume that at least one frame contains a very wideband signal
  
  tmp_decision_1 = ( tmp_lpa > ws.block.far_threshold );
  
  % -> problem: some speech may be missing.
  
  
  % ( 2 ) Consider the least conservative threshold (FRR = FRR_T)
  % Idea: in each utterance of speech, usually begin and end are weak.
  
  tmp_decision_2 = ( tmp_lpa > ws.block.frr_threshold );
  
  % -> problem: too much space to search and unwanted non-speech sources.
  % ( beamer, laptops etc. )
  
  
  % ( 3 ) Combine the two, by restricting "not-so-sure" speech detected in "tmp_decision_2"
  % to be close in space-time to "sure" speech detected in "tmp_decision_1"
  
  % Over space: sectors are assumed circular
  
  tmp_mask = tmp_decision_1;
  
  for a = 1:ws.PAR.DETECTION.SPACE_TIME_DILATION( 1 )
    
    tmp_mask = tmp_mask | tmp_decision_1( [ (a+1):end  1:a ], : );
    
    tmp_mask = tmp_mask | tmp_decision_1( [ (end-a+1):end  1:(end-a) ], : );
    
  end
  
  % Over time: we assume that time is not circular
  % which seems reasonable (at least locally! See the works of Gödel.)
  %
  % This corresponds to a buffer of "ws.PAR.SPACE_TIME_DILATION( 2 )" time frames.
  
  tmp_mask = dilation_rows( tmp_mask, ws.PAR.DETECTION.SPACE_TIME_DILATION( 2 ) );
  
  % Take the final decision
  
  n = size( ws.previous_block_ssm_log_post_act, 2 );
  
  ws.ssm_detection_decision( :, ws.block.start_frame:ws.block.stop_frame ) = ( tmp_mask( :, n+1:end ) & tmp_decision_2( :, n+1:end ) );
  
  
  % ( 4 ) Store the last few frames of this block, for the next one
  
  n = min( ws.block.stop_frame - ws.block.start_frame + 1, 3 * ws.PAR.DETECTION.SPACE_TIME_DILATION( 2 ) );
  
  ws.previous_block_ssm_log_post_act = ws.block.ssm_log_post_act( :, end-n+1:end );
  


% ----------------------------------------------------------------------
% Prepare data used by "one_array_fasttde_dirloc()"  (below)
% -> (1) where needed, compute time-domain, upsampled GCCPHAT
% -> (2) estimate time-delays

function compute_timedomain_gccphat

  global ws

  % Back to time-domain 
  % "fftshift" so that time-delay zero corresponds
  % to row "rowzero" of the "gccphat" matrix
  %
  % gccphat: 3D matrix
  % size( gccphat ) = [ 2 * ws.a_struct.framelength,  ws.a_struct.fasttde.n_pair,  ws.block.n_frames ]

  % This is only needed on frames with at least one active sector

  ind = ws.block.start_frame:ws.block.stop_frame;

  active_block_frame_ind = find( any( ws.ssm_detection_decision( :, ind ), 1 ) );

  gccphat = fftshift( real( ifft( ws.block.ssm_observed_gccphat( :, ws.a_struct.fasttde.i_pair_list, active_block_frame_ind ), 2 * ws.a_struct.framelength, 1 ) ), 1 ); 
  
  % Upsample each column of the matrix
  % Note that "up_gccphat" is a 2D matrix
  up_gccphat = upsample( gccphat( ws.a_struct.fasttde.rowstart:ws.a_struct.fasttde.rowend, : ), ws.a_struct.fasttde.upfactor );
  up_gccphat = filter( ws.a_struct.fasttde.B, 1, up_gccphat, [], 1 );
  
  % Move back to 3D   [ up_nrows,  fasttde.n_pair,  block.n_frames ]
  
  n1 = ws.a_struct.fasttde.up_nrows;
  n2 = ws.a_struct.fasttde.n_pair;
  
  ws.block.up_gccphat = zeros( [ n1  n2  ws.block.n_frames ] );
  ws.block.up_gccphat( :, :, active_block_frame_ind ) = reshape( up_gccphat, n1, n2, [] );

% ----------------------------------------------------------------------
% Localization within a sector, based on GCC-PHAT Time Delay Estimation
% (Knapp 76, Brandstein's thesis 95)
%
% Note that for each microphone pair, 
% the search for the time-delay (local maximum of GCCPHAT)
% is limited to an interval specific to this sector.
%
% Each (sector, pair)-specific interval is defined in "get_details()".

% The localization may fail, in which case "sub_az" and "sub_el" are both [].
% In case of success, both "sub_az" and "sub_el" are scalars (angles in radians).

function [ sub_az, sub_el ] = one_array_fasttde_dirloc( block_frame_ind, sector_id )

  if nargin < 2
    error( [ mfilename '.one_array_fasttde_dirloc() needs 2 input argument.' ] );
  end
  
  global ws;
  
  s = ws.ssm.sector( sector_id );

  % For each pair, limit the TDOA search range to the physically possible
  sector_up_gccphat = ws.block.up_gccphat( :, :, block_frame_ind ) + ws.a_struct.fasttde.up_range_def( :, :, sector_id );
  
  % Extract TDOA for each pair
  [ pair_gccphat_value, max_ind ] = max( sector_up_gccphat, [], 1 );

  pair_tdoa = ws.a_struct.fasttde.up_row_td( max_ind );
  
  % From this, for each subarray, deduce a direction vector,
  % i.e. azimuth and elevation in radians
  
  [ sub_az, sub_el ] = tdoa_2_subarray_direction( pair_tdoa );

  % Accept only directions within this sector
  
  sub_az = mod( sub_az, 2*pi );
  sub_el = -pi + mod( +pi + sub_el, 2*pi );
  
  ind = find( angle_within( sub_az, s.az_min, s.az_max ) & angle_within( sub_el, s.el_min, s.el_max ) );
	    
  sub_az = sub_az( ind );
  sub_el = sub_el( ind );
  
  % If more than one direction, then merge them into one
  if numel( sub_az ) > 1
    
    sub_az = angle( mean( exp( j * sub_az ) ) );
    sub_el = angle( mean( exp( j * sub_el ) ) );
        
  end
  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implementation of the "cosine" approach
% described in Branstein's thesis:
% From estimated TDOA values, for each subarray, deduce a direction vector,
% i.e. azimuth and elevation in radians

function [ sub_az, sub_el ] = tdoa_2_subarray_direction( pair_tdoa )
  
  if nargin < 1
    error( [ mfilename '.tdoa_2_subarray_direction() needs 1 input argument.' ] );
  end

  global ws

  nsubarray = ws.a_struct.fasttde.nsubarray;
  
  sub_az = repmat( NaN, nsubarray, 1 );
  sub_el = repmat( NaN, nsubarray, 1 );
  
  for j = 1:nsubarray
    
    i_i_pair_list = ws.a_struct.fasttde.i_i_pair_list( j, : );
    Pmat        = ws.a_struct.fasttde.square_subarray_Pmat_list{ j };
    
    % ( 1 ) a direction vector in the local coordinate system
    % ( the 3 cosines in Brandstein's thesis Section 7.2 )
    
    cos_alpha_j = pair_tdoa( i_i_pair_list( 1 ) ) / ws.a_struct.fasttde.maxdelay( i_i_pair_list( 1 ) );
    cos_beta_j  = pair_tdoa( i_i_pair_list( 2 ) ) / ws.a_struct.fasttde.maxdelay( i_i_pair_list( 2 ) );
    cos_gamma_j = sqrt( 1 - cos_alpha_j ^ 2 - cos_beta_j ^ 2 );
      
    if ~isreal( cos_gamma_j )
      
      % Can't apply the theory! Leave the location as "not found" -> NaNs
      sub_az( j ) = NaN;
      sub_el( j ) = NaN;
      
    else
      
      % Can apply the theory -> estimate az and el
      
      % Direction vector in the local coordinate system, in homogenous coordinates (4th row is "1")
      
      local_dir_j = [ cos_alpha_j; cos_beta_j; cos_gamma_j; 1 ];
      
      % ( 2 ) translate the direction into the global coordinate system
      
      global_dir_j = full( Pmat * local_dir_j );
      
      % Sanity check
      if abs( global_dir_j( 4 ) - 1 ) > 100 * eps
	error( [ mfilename '.tdoa_2_subarray_direction() has too many problems. Gone to the psychiatrist.' ] );
      end
      
      % ( 3 ) store azimuth and elevation (radians).
      [ sub_az( j ), sub_el( j ) ] = cart2sph( global_dir_j( 1 ), global_dir_j( 2 ), global_dir_j( 3 ) );
      
    end
    
  end
  

  
% ----------------------------------------------------------------------  
% Get the short-term past and future (azimuth location estimtes)
% This modifies: "ws.stc.past" and "ws.stc.future"
%
% Both ws.stc.past and ws.stc.future have the same format: 
% 4 by H matrix:
% - row 1:  time stamp t ( frame index )
% - row 2:  column index in "sphloc_mat"
% - row 3:  azimuth in radians
% - row 4:  cluster tag

function get_stc_past_future
  
  global ws   sphloc_mat
  
  ws.stc.past   = zeros( 4, 0 );
  ws.stc.future = zeros( 4, 0 );

  % ( 1 ) Extract past

  if ~isempty( ws.stc.last_loc_est)

    t   = sphloc_mat( 1, ws.stc.last_loc_est );
    
    ind = find( sphloc_mat( 1, 1:ws.stc.last_loc_est ) >= t - ws.PAR.STC.N_PAST_FRAMES );
    
    if ~isempty( ind )

      ws.stc.past = [ sphloc_mat( 1, ind ); ind; sphloc_mat( [ 3 7 ], ind ) ];

    end

  end

  % ( 2 ) Extract future

  if ~isempty( ws.stc.last_loc_est )
    
    ind = ws.stc.last_loc_est+1:size( sphloc_mat, 2 );

  else

    ind = 1:size( sphloc_mat, 2 );
    
  end
  
  % Limit the complexity
  % (see IDIAP RR 04-14)

  ind = ind( 1:min( numel( ind ), ws.PAR.STC.N_FUTURE_SAMPLES ) );

  if ~isempty( ind )

    ws.stc.future = [ sphloc_mat( 1, ind ); ind; sphloc_mat( [ 3 7 ], ind ) ];

  end
  
    
% ----------------------------------------------------------------------
% Apply short-term clustering
% This modifies "ws.stc.future" (row 4: cluster tags that are NaNs become non-NaN integers)
% and "ws.stc.last_tag"

function short_term_clustering

  global ws

  np = size( ws.stc.past, 2 );
  nf = size( ws.stc.future, 2 );
  n = np + nf;

  % Anything to do?
  if nf < 1
    % No
    return;
  end

  % Yes

  if nf == 1

    %%%%%%%%%%
    % Special case: if there is only one sample to add,
    % we can do things very simply (only step 3 is needed,
    % as defined in IDIAP RR 04-14).
    %
    % This special case was implemented purely for speed purposes.
    
    delta_t  = abs( ws.stc.past( 1, : ) - ws.stc.future( 1, 1 ) );

    ind = find( delta_t <= ws.PAR.STC.T_SHORT );

    % Are we sure that we need to create a new cluster?
    if isempty( ind )
      % Yes
      % Create a new cluster
      ws.stc.last_tag       = ws.stc.last_tag + 1;
      ws.stc.future( 4, 1 ) = ws.stc.last_tag;
      % Done
      return;
    end

    % No
    % Test several hypotheses

    delta_t  = delta_t( ind );
    delta_az = -pi + mod( ws.stc.past( 3, ind ) - ws.stc.future( 3, 1 ) + pi, 2 * pi );
    past_tag = ws.stc.past( 4, ind );
    
    past_tag_list  = unique( past_tag );
    n_past_cluster = numel( past_tag_list );
    
    % Evaluate pair-wise likelihood involving the new sample
    
    past_log_lkld_same    = log_normpdf( delta_az, 0, ws.stc.model.same( delta_t + 1 ) );
    past_log_lkld_notsure = log_normpdf( delta_az, 0, ws.stc.model.notsure( delta_t + 1 ) );
    
    % All we need to do is to evaluate (n_past_cluster + 1) hypotheses
    % using the criterion defined in IDIAP RR 04-14.
    % The "+ 1" below corresponds to creating a new cluster.
    %
    % In fact we only evaluate the terms in the criterion
    % that involve the new sample, since other terms remain unchanged.
    
    hyp_log_lkld = repmat( NaN, 1, n_past_cluster + 1 );
    
    for a = 1:n_past_cluster

      % Estimate the likelihood of the data when putting the new sample
      % into one of the existing past clusters
      
      ind_same    = find( past_tag == past_tag_list( a ) );
      ind_notsure = setdiff( 1:numel( past_tag ), ind_same );

      hyp_log_lkld( a ) = sum( past_log_lkld_same( ind_same ) ) + sum( past_log_lkld_notsure( ind_notsure ) );
      
    end
    
    % Estimate the likelihood of the data when creating a new cluster,
    % with just the new sample in it.
    
    hyp_log_lkld( n_past_cluster + 1 ) = sum( past_log_lkld_notsure );

    % Sanity check

    if any( isnan( hyp_log_lkld ) )
      error( [ mfilename '.short_term_clustering() is NaN-insane.' ] );
    end

    % Now select the most likely hypothesis

    [ tmp, max_ind ] = max( hyp_log_lkld );

    if max_ind <= n_past_cluster

      % Put the new sample into an existing cluster
      ws.stc.future( 4, 1 ) = past_tag_list( max_ind );
      
    else
    
      % Create a new cluster
      ws.stc.last_tag       = ws.stc.last_tag + 1;
      ws.stc.future( 4, 1 ) = ws.stc.last_tag;
      
    end

    % Done with the special case
    return;

  end

  %%%%%%%%%%  
  % General case: more than one sample to add.
  % Both steps 2. and 3. are fully implemented.
  % (see IDIAP RR 04-14).

  % ( 0 ) Get data
  
  past_and_future = [ ws.stc.past  ws.stc.future ];
  
  % Frame index
  % we don't need sphloc_mat column index here
  t   = past_and_future( 1,: );

  az  = past_and_future( 3,: );  % azimuth in radians
  tag = past_and_future( 4,: );  % Cluster tag (an integer, future samples have NaNs)
  
  % ( 1 ) Get the graph corresponding to past samples
  
  % "full_graph" will eventually contain both past and future
  % "1" indicate an edge between two vertices (row) and (col)
  %
  % Only the strictly superior triangular part is used.
  
  full_graph = sparse( n, n );
  
  full_graph( 1:np, 1:np ) = triu( repmat( tag( 1:np ), np, 1 ) == repmat( tag( 1:np ).', 1, np ), 1 );
  
  
  % ( 2 ) Get the "same" likelihood and "different" likelihood
  
  delta_t  = abs( repmat( t, n, 1 ) - repmat( t(:), 1, n ) );
  delta_az = -pi + mod( pi + repmat( az, n, 1 ) - repmat( az(:), 1, n ), 2 * pi );
  
  old_warning_state = warning; warning off;
  
  ind = find( delta_t <= ws.PAR.STC.T_SHORT );
  
  tmp_x     = [ delta_az( ind ) ];  % Observed data:   numel( ind ) by 1 matrix

  % "same object"
  % put NaN values for "undefined" likelihoods,
  % i.e. whenever the delta_t is larger than ws.PAR.STC.T_SHORT.
  
  loglike_same = repmat( NaN, n, n );
  
  loglike_same( ind ) = log_normpdf( tmp_x, 0, ws.stc.model.same( delta_t( ind ) + 1 ).' );
  
  % "not sure"
  % keep NaN values for "undefined" likelihoods,
  % i.e. whenever the delta_t is larger than ws.PAR.STC.T_SHORT.
  
  loglike_notsure = repmat( NaN, n, n );
  
  loglike_notsure( ind ) = log_normpdf( tmp_x, 0, ws.stc.model.notsure( delta_t( ind ) + 1 ).');
  
  warning( old_warning_state );
  

  % ( 3 ) On future frames only,
  % get the likelihoods 
  % get a list of ALL possible graphs
  
  future.loglike_notsure = loglike_notsure( 1+np:end, 1+np:end );
  future.loglike_same    = loglike_same( 1+np:end, 1+np:end );
  
  future.graph_list      = ws.stc.all_graphs( nf ).graph_list;
  
  % Pick the best one
  % ( Step 2. in IDIAP RR 04-14 )
  
  [ best_ind, graph_list ] = select_graph( future.loglike_same, future.loglike_notsure, future.graph_list );
  future.graph_list = graph_list;
  
  new_graph = future.graph_list{ best_ind };
  
  siz = size( new_graph );
  
  % Sanity check
  if ( numel( siz ) ~= 2 ) | ( siz(1) ~= siz(2) )
    error( [ mfilename '.short_term_clustering(): check for sanity!' ] );
  end
  
  % Sanity check
  if siz(1) ~= nf
    error( [ mfilename '.short_term_clustering(): check for sanity again!' ] );
  end
  
  full_graph( 1+np:end, 1+np:end ) = new_graph;
  

  % ( 4 ) Connect past and future
  % ( Step 3. in IDIAP RR 04-14 )

  if np > 0

    % Construct all possible merges 
    % note that the list includes the "no merge" case
    
    r = 1:np;
    c = np+1:n;
    
    graph_1 = full_graph( r, r );
    graph_2 = full_graph( c, c );
    
    % Limitation of the permitted "same" edges
    % permitted_1 is a N x N matrix of binary values,
    % where "1" on an edge indicates a permitted "same" decision for this edge,
    % and   "0" on an edge indicates a forbidden "same" decision for this edge,
    % 
    % -> prune when log lkld ratio <= ws.PAR.STC.LLR_THRESHOLD
    permitted_1 = double( ( loglike_same( r, c ) - loglike_notsure( r, c ) ) > ws.PAR.STC.LLR_THRESHOLD );
    permitted_1( find( isnan( loglike_same( r, c ) ) ) ) = NaN;
    
    graph_list = construct_merge_two( graph_1, graph_2, permitted_1, ws.PAR.STC.N_MERGES_MAX );
    
    % Pick the best one in a ML sense
    
    [ best_ind, graph_list ] = select_graph( loglike_same, loglike_notsure, graph_list );	
    
    % Store the merge
    
    full_graph( r, c ) = graph_list{ best_ind }( r, c );
    
    % Sanity check
    if any( full_graph(:) ~= graph_list{ best_ind }(:) )
      
      error( [ mfilename '.short_term_clustering(): I do not understand what is happening.' ] );
      
    end
    
  end
  
  
  % ( 5 ) Return the result: convert the graph back to a list of tags
  
  for a = 1:n
    
    % Make sure this sample has a tag
    if isnan( tag( a ) )
      
      ws.stc.last_tag = ws.stc.last_tag + 1;
      tag( a ) = ws.stc.last_tag;
      
    end

    % Propagate the tag to other samples in the same cluster
    tag( find( full_graph( a, : ) ) ) = tag( a );
    
    % Are we done?
    if all( ~isnan( tag ) )
      break;
    end
    
  end
  
  % Sanity check

  if any( ws.stc.past( end, : ) ~= tag( 1:np ) )
    error( [ mfilename '.short_term_clustering() is not supposed to change tags of past frames.' ] );
  end
  
  % Finally update the tags in the workspace

  ws.stc.future( end, : ) = tag( np+1:end );

  
% ----------------------------------------------------------------------
% Tool used by the "short_term_clustering()" function
% Select the ML graph out of list,
% given "same" likelihoods and "different" likelihoods.
% (IDIAP RR 04-14)

function [ best_ind, graph_list, graph_like, best_h1_post ] = select_graph( loglike_same, loglike_notsure, graph_list )

  n = size( loglike_same, 1 );
  ng = length( graph_list );
  
  % ( 1 ) Compute graph likelihoods
  
  graph_like = repmat( NaN, 1, ng );
  
  % Likelihood computation is restricted to |t_i - t_j| <= T_{short}
  bin_mask = ~isnan( loglike_same );
  
  for a = 1:ng

    g = graph_list{ a };
    graph_like( a ) = sum( loglike_same( find( g & bin_mask ) ) ) + sum( loglike_notsure( find( ~g & bin_mask ) ) );
    
  end

  % Pick the ML one
  [tmp, best_ind] = max( graph_like );

  % Do we need posteriors?

  if nargout < 4

    % No
    return;

  end

  % Yes

  % ( 2 ) Compute posteriors for each possible local decision,
  % using all possible graphs.
  %
  % This can be useful for confident clustering (IDIAP RR 04-14)
  % ...but definitely useless in audio meeting recordings (e.g. M4 corpus).
  
  graph_scaled = graph_like - logsum( graph_like );
  
  h0_post = repmat( NaN, n, n );
  h0_post( find( bin_mask ) ) = -Inf;
  
  h1_post = repmat( NaN, n, n );
  h1_post( find( bin_mask ) ) = -Inf;
  
  for a = 1:ng
  
    g = graph_list{ a };
    
    ind = find( ~g & bin_mask );
    ind = ind(:).';
    m = h0_post( ind );
    m( 2,: ) = graph_scaled( a );
    h0_post( ind ) = logsum( m );
    
    ind = find( g & bin_mask );
    ind = ind(:).';
    m = h1_post( ind );
    m( 2,: ) = graph_scaled( a );
    h1_post( ind ) = logsum( m );
    
  end
  
  best_h1_post = h1_post;
  best_h1_post( find( ~graph_list{ best_ind } ) ) = 0;
  best_h1_post = sparse( best_h1_post );
  
  

% ----------------------------------------------------------------------
% Tool used by the "short_term_clustering()" function
% to construct all possible merges of two graphs.
%
% ( used to select the best merge of the past window and the future window ).
%
% in_graph1 and in_graph2 are N x N matrices of binary values.
%
% Optional limitation of the permitted "same" edges:
% permitted_1 is a N x N matrix of binary values,
% where "1" on an edge indicates a permitted NEW "same" decision for this edge,
% and   "0" on an edge indicates a forbidden NEW "same" decision for this edge.
%
% n_merges_max limits the number of constructed graphs (default: +Inf).
% 


function graph_list = construct_merge_two( in_graph1, in_graph2, permitted_1, n_merges_max, n_elements_min )
  
  % Deal with parameters
  
  if nargin < 3
    error( [ mfilename 'construct_merge_two() needs at least 3 input arguments' ] )
  end
  
  if ~exist( 'n_merges_max', 'var' )
    n_merges_max = +Inf;
  end
 
  if ~exist( 'n_elements_min', 'var' )
    n_elements_min = 1;
  end
  
  % Quick check
  siz1 = size( in_graph1 );
  siz2 = size( in_graph2 );
  siz3 = size( permitted_1 );
  
  if siz1(1)~=siz1(2)
    error( [ mfilename 'construct_merge_two(): "in_graph1" must be a square matrix !' ] );
  end
  n1 = siz1( 1 );
  
  if siz2(1)~=siz2(2)
    error( [ mfilename 'construct_merge_two(): "in_graph2" must be a square matrix !' ] );
  end
  n2 = siz2( 1 );
  
  if any( [ n1 n2 ] ~= siz3 )
    error( [ mfilename 'construct_merge_two(): wrong size for "permitted_1" !' ] );
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % If nothing to do, just return
  
  if ~any( permitted_1(:) == 1 )
    
    new_graph = in_graph1;
    n = 1:n2;
    new_graph( end+n, end+n ) = in_graph2;
    
    graph_list{ 1 } = new_graph;
    
    return;
  
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Otherwise get to work
  % First, convert graphs to tags
  
  tag1 = graph_to_tag( in_graph1 );
  tag2 = max( tag1 ) + graph_to_tag( in_graph2 );
  
  tag1_list = unique( tag1 );
  tag2_list = unique( tag2 );

  ntag1 = length( tag1_list );
  ntag2 = length( tag2_list );

  % List the clusters we refuse to merge with any other cluster
  % Elements of "exclusion1_list" are indices in "tag1_list"
  % Elements of "exclusion2_list" are indices in "tag2_list"
  exclusion1_list = [];
  exclusion2_list = [];
  
  if n_elements_min > 1
    
    for a = 1:length( tag1_list )
      a_tag = tag1_list( a );
      n_elements = sum( tag1 == a_tag );

      if n_elements < n_elements_min
	exclusion1_list( end+1 ) = a;
      end
    
    end
    
    for a = 1:length( tag2_list )
      a_tag = tag2_list( a );
      n_elements = sum( tag2 == a_tag );

      if n_elements < n_elements_min
	exclusion2_list( end+1 ) = a;
      end
    
    end
    
  end
  
  % Convert permitted_1 into a list of excluded merges
  %
  % forbidden_merge is a cell array of ntag2 vectors of indices,
  % each index is an index of an element the tag1 vector.
  
  forbidden_merge = cell( 1, ntag2 );
  
  m = zeros( ntag1, ntag2 );
  
  [u,v] = find( permitted_1 == 0 );
  
  if ~isempty( u );
    
    % Convert to tag 
    u = tag1( u );
    v = tag2( v );
    
    % Convert to index in each tag list
    tag1_to_index1( tag1_list ) = 1:ntag1;
    tag2_to_index2( tag2_list ) = 1:ntag2;
    
    u = tag1_to_index1( u );
    v = tag2_to_index2( v );
    
    % Mark the forbidden merges
    m( sub2ind( size( m ), u, v ) ) = 1;
    
  end
  
  for a = 1:ntag2
   
    % For each cluster in graph2
    % we list the forbidden clusters in graph1
    
    forbidden_merge{ a } = find( m( :, a ) );
    
  end
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Second, define all possible merges between clusters
  % of in_graph1 and clusters of in_graph2. 
  % ( No merge between clusters of the same graph. )
  % 
  % Recursive implementation:
  % Each element of "merge_list" is a ntag1 x ntag2 sparse matrix of zeroes and ones,
  % such that there is only one "one" on each row.
  % We recursively build a list of (n+1)-column matrices from each n-column matrix.
  % 
  % I have no time for a better explanation
  % If this is not clear: lathoud@idiap.ch
  
  merge_list = { sparse( [], [], [], ntag1, 0 ) };

  reached_max = 0;
  
  while ~reached_max  &  size( merge_list{ 1 }, 2 ) < ntag2     % Do we have the right number of columns?
    
    % If not, then add a column : store all possibilities
    % that respect the constraint of only one "one" on each row
    
    new_merge_list = {};
    
    for a = 1:length( merge_list )
      
      if length( new_merge_list ) >= n_merges_max
	s = warning;
	if ~strcmp( s( 1 ).state, 'off' )
	 disp( 'construct_merge_two: WARNING! Reached the maximum number of merges.' ); 
	end
	
	% Fall back on the existing list of merges
	reached_max = 1;
	break;
	
      end
      
      m = merge_list{ a };
      
      % We can always add a column of zeroes (means no additional merge)
      
      new_merge_list{ end+1 } = [ m  zeros( ntag1, 1 ) ];
      
      if ~ismember( size( m, 2 ) + 1, exclusion2_list )
      
	% All possible non-zero columns
	
	if size( m, 2 ) < 1
	  ind = 1:size( m, 1 );
	else
	  ind = find( ~any( m, 2 ) );
	end
	
	ind = setdiff( ind, [ exclusion1_list  forbidden_merge{ size( m, 2 ) + 1 } ] );
	
	for b = ind(:).'
	  new_merge_list{ end+1 } = [ m  (1:ntag1).'== b ];
	end
	
      end
      
    end

    if reached_max
      break;
    end
    
    % Update with the new list of matrices with n+1 columns
    
    merge_list = new_merge_list;
    
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Third convert back to graphs
  % Each vector of tags becomes a matrix of binary values.

  ind = [];
  
  graph_list = {};
  
  for a = 1:length( merge_list )
    
    [u,v] = find( merge_list{ a } );
    v = tag2_list( v );
    
    m = ( repmat( tag2(:), 1, length( v ) )  ==  repmat( v(:).', length( tag2 ), 1 ) );
    [ r, c ] = find( m );
    
    tag2bis = tag2;
    tag2bis( r ) = u( c );
    
    tag = [ tag1 tag2bis ];
    
    graph_list{ a } = tag_to_graph( tag );
    
    if isempty( ind )
      m = zeros( size( graph_list{ a } ) );
      m( 1:n1, n1+1:end ) = permitted_1;
      ind = find( isnan( m ) );
    end
    
    graph_list{ a }( ind ) = 0;
    
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Fourth make sure the constructed graphs
  % don't have the "extra ones"
  % ( e.g. for T > T_short )
  
  if ~all( permitted_1 == 1 )
    
    for a = 1:length( graph_list )
      graph_list{ a }(1:n1, n1+1:end ) = graph_list{ a }(1:n1, n1+1:end ) & ( permitted_1 == 1 );
    end
    
  end


% --------------------------------------------------
% Convert a N x N graph (matrix of binary values)
% into a vector of N tag values.

function tag = graph_to_tag( graph )
  
  if nargin < 1
    error( [ mfilename '.graph_to_tag() needs one input argument' ] );
  end
  
  N = size( graph, 1 );
  n_tags = 0;
  
  tag = zeros( 1, N );
  
  for a = 1:N
    
    % Make sure this vertice has a tag
    if ~tag( a )
      n_tags = n_tags + 1;
      tag( a ) = n_tags;
    end
    
    % Propagate this tag to connected vertices
    tag( find( graph( a, : ) ) ) = tag( a );
    
  end

% --------------------------------------------------
% Convert a vector of N tag values
% into a N x N graph (matrix of binary values)

function graph = tag_to_graph( tag )
  
  if nargin < 1
    error( [ mfilename '.tag_to_graph() needs one input argument' ] );
  end
  
  n = length( tag );
  graph = sparse( triu( repmat( tag(:), 1, n ) == repmat( tag(:).', n, 1 ), 1 ) );
  

  
% ----------------------------------------------------------------------
% Store the result of short-term clustering
% Input:  "ws.stc.future"
% Output: this modifies: "sphloc_mat", "ws.stc.last_loc_est"

function update_stc_past_future
  
  global ws   sphloc_mat
  
  % ( 1 ) Update "sphloc_mat"
  % -> we only have to replace the NaNs with cluster tags

  ind = ws.stc.future( 2,: );

  if any( ~isnan( sphloc_mat( 7, ind ) ) )
    error( [ mfilename '.update_stc_past_future() is non-NaN-insane.' ] );
  end
  
  sphloc_mat( 7, ind ) = ws.stc.future( 4, : );
  
  % ( 2 ) Update various counters in "ws.stc"
  
  ws.stc.last_frame   = max( ws.stc.future( 1, : ) );
  
  ws.stc.last_loc_est = max( ws.stc.future( 2, : ) );

  ws.stc.n_new_sample = ws.stc.n_new_sample - size( ws.stc.future, 2 );


% ----------------------------------------------------------------------
% Reduce the amount of SAM-SPARSE-MEAN data
% Input:  "ws.block.previous_ssm_activeness"
%         "ws.block.ssm_activeness"
% Output: "ws.block.reduced_ssm_activeness"
% ( useful for EM )

function reduce_ssm_data
  
  global ws

  % Note: the code is replicated twice,
  % in order to avoid data replication.
  
  % ( 1 ) Reduce "previous data"
  
  n_frames     = size( ws.block.previous_ssm_activeness, 2 );
  max_n_frames = ws.a_struct.datarepr_maxprevious_frames; 
  
  if n_frames <= max_n_frames
    
    % Nothing to do
    
    ws.block.reduced_ssm_activeness = ws.block.previous_ssm_activeness;
    
  else
    
    % Reduction is necessary
    
    [ tmp, ind ] = sort( max( ws.block.previous_ssm_activeness, [], 1 ) );
    
    ind = ind( min( n_frames, 1 + round( ( n_frames-1 ) * ( 0.5:max_n_frames ) / max_n_frames ) ) );
    
    ws.block.reduced_ssm_activeness = ws.block.previous_ssm_activeness( :, ind );
    
  end
  
  
  % ( 2 ) Reduce "current data"
  
  n_frames     = size( ws.block.ssm_activeness, 2 );
  max_n_frames = ws.a_struct.datarepr_maxcurrent_frames; 
  
  if n_frames <= max_n_frames
    
    % Nothing to do
    
    ws.block.reduced_ssm_activeness = [ ws.block.reduced_ssm_activeness  ws.block.ssm_activeness ];
    
  else
    
    % Reduction is necessary
    
    [ tmp, ind ] = sort( max( ws.block.ssm_activeness, [], 1 ) );
    
    ind = ind( min( n_frames, 1 + round( ( n_frames-1 ) * ( 0.5:max_n_frames ) / max_n_frames ) ) );
    
    ws.block.reduced_ssm_activeness = [ ws.block.reduced_ssm_activeness  ws.block.ssm_activeness( :,ind ) ];
    
  end
  
  % Sanity check
  if size( ws.block.reduced_ssm_activeness, 2 ) > ws.a_struct.datarepr_maxprevious_frames + ws.a_struct.datarepr_maxcurrent_frames
    error( [ mfilename '.reduce_ssm_data(): insanity detected.' ] );
  end
  
  
% ----------------------------------------------------------------------
% Unsupervised fit of a detection model on "ws.block.reduced_ssm_activeness".
% The chosen model is 2 Dirac + 2 Gammas + 1 Shifted Erlang.
% See IDIAP RR 05-52 for theory and implementation details.
%
% Note that we use a Shifted Erlang instead of a Shifted Rice,
% because it is more robust for smaller amounts of data.
%
% For a detailed description of what each parameter name means,
% see "compute_post_2d2g1se.m".

function fit_model
  
  global ws 
  
  % This is not supposed to be changed
  GAMMA0_FLOOR = 1 + 1e-10;
 
  
  ws.z_ind  = find( ws.block.reduced_ssm_activeness == 0 );
  ws.nz_ind = find( ws.block.reduced_ssm_activeness ~= 0 );
 
  ws.log_data_nz_ind = log( ws.block.reduced_ssm_activeness( ws.nz_ind ) );
  
  % ( 1 ) Automatic initialization
  % We test several initialization, where for each one
  % the data is splitted in two (above and below a threshold)
      
  % Flat (1-D), ordered representation of the data

  sort_data = sort( ws.block.reduced_ssm_activeness(:) );
  
  % List of thresholds that we are going to try
  
  a_min = sort_data( 1 );
  a_max = sort_data( end );
  
  if ws.PAR.DETECTION.AUTO_INIT_N < 2
  
    % Only one threshold (I don't advise this)

    thr_list = ( a_max + a_min ) / 2;
    
  else
    
    % Many thresholds
    
    n1 = max( 1, round( ws.PAR.DETECTION.AUTO_INIT_N / 2 ) );
    
    n2 = max( 0, ws.PAR.DETECTION.AUTO_INIT_N - n1 );
    
    % Multiple initializations
    % (we'll pick the ML one)
    % -> this way we avoid gross mistakes
    %    that would lead to bad local optima
        
    thr_list = a_min + (0.5:n1) / n1 * ( a_max - a_min );
    
    % Add the same number of thresholds,
    % regularly spaced along percentiles
    
    n = numel( sort_data );
    ind = 1 + round( ( 0.5:n2 ) / n2 * ( n - 1 ) );
    
    tmp = sort_data( ind );
    
    thr_list = [ thr_list(:).' tmp(:).' ];
     
  end

  thr_list = sort( thr_list );
  
  for a = 1:numel( thr_list )
    
    thr = thr_list( a );
    
    % Gamma pdfs: moment method, same parameters for both
    
    ind = find( sort_data <= thr );
    
    % We need at least two samples to apply the moment method
    a_mean = mean( sort_data( ind ) );
    a_var  = var( sort_data( ind ) );
	
    if ( numel( ind ) < 2 ) | (  a_var < eps )
      
      log_plkld_list( a ) = -Inf;
      
      a_model.a0   = [];
      a_model.b0   = [];
      a_model.w00d = [];
      a_model.w01d = [];
      a_model.gamma00 = [];
      a_model.beta00  = [];
      a_model.gamma01 = [];
      a_model.beta01  = [];
      a_model.beta11 = [];
      
    else
      
      % Ensure that GammaPdf( 0, gamma00, beta00 ) = 0
      
      a_model.gamma00 = max( GAMMA0_FLOOR, a_mean ^ 2 / a_var );
      a_model.beta00  = max( 100*eps, a_mean / a_model.gamma00 );
      
      % Same parameters for ( gamma01, beta01 )
      
      a_model.gamma01 = a_model.gamma00;
      a_model.beta01  = a_model.beta00;
      
      % Shifted Erlang pdf: moment method as well
      
      t11 = a_model.gamma01 * a_model.beta01;
      grt_thr = max( t11, thr );
      grt_ind = find( sort_data > grt_thr );
      if isempty( grt_ind )
	error( [ mfilename ' could not find data to initialize the shifted Rice!' ] );
      end
      
      % Implicitely: gamma11 = 2
      a_model.beta11 = mean( sort_data( grt_ind ) - t11 ) / 2;
      
      % We are almost sure to get those crappy warnings
      % so let's shut if off during initilization
      old_warning_state = warning;
      warning off;
      
      
      % Update priors
      
      tmp_max = max( ws.block.reduced_ssm_activeness, [], 1 );
      
      tmp_max_thr = exp(mean(log(tmp_max(find(tmp_max)))));
      
      a_model.a0 = min( 0.9, max( 1e-1, mean( tmp_max <= tmp_max_thr ) ) );
            
      
      tmp = ws.block.reduced_ssm_activeness( :, tmp_max > tmp_max_thr );
      
      a_model.b0 = min( 0.9, max( 0.1, mean( tmp(:) <= grt_thr ) ) ); 
      
      
      % Update mixture weights
      
      tmp = ws.block.reduced_ssm_activeness( :, tmp_max <= tmp_max_thr );
      
      a_model.w00d = min( 0.9, max( ws.PAR.DETECTION.WD_WG_FLOOR, mean( tmp(:) == 0 ) ) );
     
      
      tmp = ws.block.reduced_ssm_activeness( :, tmp_max > tmp_max_thr );
      
      a_model.w01d = min( 0.9, max( ws.PAR.DETECTION.WD_WG_FLOOR, mean( tmp(:) == 0 ) ) );

      
      % Compute lkld
      
      try	
	[ tmp, tmp, log_plkld_list( a ) ] = compute_post_2d2g1se( a_model, ws.block.reduced_ssm_activeness );
      catch
	% compute_post_2d2g1se failed (not supposed to happen)
	log_plkld_list( a ) = -Inf;
      end
      
      warning( old_warning_state );
      
    end
    
    model_list( a ) = a_model;
    
  end  % loop through thresholds
     
  % Pick the ML
  
  [ tmp, max_ind ] = max( log_plkld_list );
  
  a_model = model_list( max_ind );
  
  if ws.PAR.DETECTION.DEBUG
    
    figure; plot(thr_list,log_plkld_list,'.-')
    xlabel( 'init threshold' ); ylabel( 'alldata log plkld' );
    hold on; plot( thr_list( max_ind ), log_plkld_list( max_ind ), 'ro' );
    title( strrep( [ mfilename ' initialization' ], '_', '\_' ) );
    keyboard
    
  end
  
  
  % ( 2 ) EM
      
  ws.best    = [];
  
  ws.iter = 0;
  ws.finished = ( ws.iter >= ws.PAR.DETECTION.MAX_ITER );
  ws.previous_alldata_log_plkld = [];
  
  a_new.model = a_model;
  
  
  while ~ws.finished
    
    if ws.PAR.DETECTION.VERBOSE > 1
      tic
    end
    
    % E step: estimate posteriors
    
    [ frame_log_post, sector_log_post, a_new.alldata_log_plkld ] = compute_post_2d2g1se( a_new.model, ws.block.reduced_ssm_activeness );
    
    % Verbosity
    
    if ws.PAR.DETECTION.VERBOSE > 1
      
      disp( ' ' );
      
      disp( sprintf( [ mfilename ': model at iteration ws.iter:%d' ], ws.iter ) );
      disp( a_new.model );
      
      previous_alldata_log_plkld = ws.previous_alldata_log_plkld;
      
      if isempty( previous_alldata_log_plkld )
	disp( sprintf( [ mfilename ': alldata_log_plkld:%g' ], a_new.alldata_log_plkld ) );
      else
	disp( sprintf( [ mfilename ': alldata_log_plkld:%g (variation:%g)' ], a_new.alldata_log_plkld, a_new.alldata_log_plkld - previous_alldata_log_plkld ) );
      end
      
      ws.previous_alldata_log_plkld = a_new.alldata_log_plkld;
      
    end
    
    
    % Look for improvement
    
    % Initialization
    is_improving = isempty( ws.best );

    if ~is_improving
      % General case: compare partial likelihoods
      is_improving = ( ws.best.alldata_log_plkld < a_new.alldata_log_plkld - ws.PAR.DETECTION.LKLD_CONV );
    end
    
    % If improvement, then update the "best so far"
    if is_improving
      ws.best = a_new;
    end
    
    
    % Termination test
    ws.finished = ( ( ~is_improving ) & ( ws.iter >= ws.PAR.DETECTION.MIN_ITER ) ) | ( ws.iter >= ws.PAR.DETECTION.MAX_ITER );
    
    % M step: update parameters

    if ~ws.finished

      % Update priors of silent frame and active frame
      
      tmp = my_logsum_fast( [ reshape( frame_log_post( :,:,1 ), [], 1 ) reshape( frame_log_post( :,:,2 ), [], 1 ) ] );
      tmp = exp( tmp - my_logsum_fast( tmp ) );
      
      a_new.model.a0 = tmp( 1 ) / sum( tmp );  % For minor numerical errors
      a_new.model.a0 = max( ws.PAR.DETECTION.A0_A1_FLOOR, min( 1 - ws.PAR.DETECTION.A0_A1_FLOOR, a_new.model.a0 ) );
      
      % Update priors of silent sector and active sector, given that the frame is active
      tmp = my_logsum_fast( [ reshape( sector_log_post( :,:,2 ), [], 1 ) reshape( sector_log_post( :,:,3 ), [], 1 ) ] );
      tmp = exp( tmp - my_logsum_fast( tmp ) );
      
      a_new.model.b0 = tmp( 1 ) / sum( tmp );  % For minor numerical errors
      a_new.model.b0 = max( ws.PAR.DETECTION.B0_B1_FLOOR, min( 1 - ws.PAR.DETECTION.B0_B1_FLOOR, a_new.model.b0 ) );
      
      % Update weights of the ( silent sector, silent frame ) mixture
      tmp_mat = sector_log_post( :,:,1 );

      % Note that often, this one is empty (absence of zero data) -> logsum is -Inf by definition
      tmp      = my_logsum_fast( reshape( tmp_mat( ws.z_ind  ), [], 1 ) );
      
      tmp( 2 ) = my_logsum_fast( reshape( tmp_mat( ws.nz_ind ), [], 1 ) );
      
      tmp = exp( tmp - my_logsum_fast( tmp ) );
      a_new.model.w00d = tmp( 1 ) / sum( tmp );  % For minor numerical errors
      a_new.model.w00d = min( 1 - ws.PAR.DETECTION.WD_WG_FLOOR, max( ws.PAR.DETECTION.WD_WG_FLOOR, a_new.model.w00d ) ); % Flooring
      
      % Update weights of the ( silent sector, active frame ) mixture
      tmp_mat = sector_log_post( :,:,2 );
      
      % Note that often, this one is empty (absence of zero data) -> logsum is -Inf by definition
      tmp      = my_logsum_fast( reshape( tmp_mat( ws.z_ind  ), [], 1 ) );
      
      tmp( 2 ) = my_logsum_fast( reshape( tmp_mat( ws.nz_ind ), [], 1 ) );
      
      tmp = exp( tmp - my_logsum_fast( tmp ) );
      a_new.model.w01d = tmp( 1 ) / sum( tmp );  % For minor numerical errors
      a_new.model.w01d = min( 1 - ws.PAR.DETECTION.WD_WG_FLOOR, max( ws.PAR.DETECTION.WD_WG_FLOOR, a_new.model.w01d ) ); % Flooring
      
      
      % Update the ( gamma00, beta00 ) parameters
      
      t00 = a_new.model.gamma00 * a_new.model.beta00;
      
      X0 = [ t00;
	     a_new.model.gamma00;
	     NaN ];

      ws.negloglike_w0 = exp( sector_log_post( :,:,1 ) );
      ws.negloglike_w1 = [];
      
      fval0 = negloglike( X0 );
      
      [ X, fval ] = fminsearch( @negloglike, X0, ws.PAR.DETECTION.FMINSEARCH_OPTIONS );
      
      % Sanity check
      if -fval < -fval0
	error( [ mfilename ' decreased the likelihood during ( sigma00, V00 ) optimization!' ] );
      end
     
      a_new.model.gamma00 = X( 2 );
      a_new.model.beta00  = X( 1 ) / X( 2 );
            
      % Update the ( gamma01, beta01, beta11 ) parameters
      
      t01 = a_new.model.gamma01 * a_new.model.beta01;

      X0 = [ t01;
	     a_new.model.gamma01;
	     a_new.model.beta11 ];
      
      ws.negloglike_w0 = exp( sector_log_post( :,:,2 ) );
      ws.negloglike_w1 = exp( sector_log_post( :,:,3 ) );
      
      
      fval0 = negloglike( X0 );
      
      [ X, fval ] = fminsearch( @negloglike, X0, ws.PAR.DETECTION.FMINSEARCH_OPTIONS );
      
      if -fval < -fval0
	error( [ mfilename ' decreased the likelihood during ( sigma01, V01, beta11 ) optimization !' ] );
      end
      
      a_new.model.gamma01 = X( 2 );
      a_new.model.beta01  = X( 1 ) / X( 2 );
      
      a_new.model.beta11  = X( 3 );
      
    end      
    
    
    % Prepare for the next iteration
    ws.iter = ws.iter + 1;
    
    if ws.PAR.DETECTION.VERBOSE > 1
      disp( ' ' );
      disp( sprintf( [ mfilename ': iter duration was %g sec' ], toc ) );
    end
    
    if ws.PAR.DETECTION.DEBUG
      disp( [ mfilename ': debug mode, "keyboard" at the end of each iteration.' ] );
      keyboard
    end
    
  end  % of "while" loop


  
  % ( 3 ) Save a bit of memory
  
  ws.z_ind = [];
  ws.nz_ind = [];
  ws.log_data_nz_ind = [];
  
  
  % ( 4 ) Return the actual result
  
  ws.block.model = ws.best.model;
  
  % Also return log posteriors of activity in a ( frame, sector )
  
  [ tmp, sector_log_post ] = compute_post_2d2g1se( ws.block.model, ws.block.ssm_activeness );
  
  ws.block.ssm_log_post_act = sector_log_post( :,:,3 );

  % As well as the log posteriors of the reduced data

  [ tmp, sector_log_post ] = compute_post_2d2g1se( ws.block.model, ws.block.reduced_ssm_activeness );

  ws.block.reduced_ssm_log_post_act = sector_log_post( :,:,3 );
  
  % Verbosity to conclude the training
  
  if ws.PAR.DETECTION.VERBOSE > 1
    disp( ' ' );
    disp( sprintf( [ mfilename ': training finished.' ] ) )
    
    disp( sprintf( [ '  ws.iter:%d, ws.best.alldata_log_plkld:%g, ws.best.model:' ], ws.iter, ws.best.alldata_log_plkld ) );
    disp( ws.best.model );

  end
  
  
  
  %%% Debug: if activated, show some plots
  
  if ws.PAR.DETECTION.DEBUG

    model = ws.best.model;
    
    % Show the distributions
    
    [h.w, h.x] = hist( ws.block.reduced_ssm_activeness(:), ws.PAR.DETECTION.DEBUG_HIST_NBINS );
    h.x = h.x(:);
    h.w = running_mean( h.w, 3 );
    h.w = h.w(:) / sum( h.w );
    
    [ h.frame_log_post, h.sector_log_post, h.alldata_log_plkld, h.frame_log_plkld, h.frame_log_pjoint, h.sector_log_pjoint_cond ] = compute_post_2d2g1se( model, h.x );
    h.sector_log_pjoint = cat( 3, h.sector_log_pjoint_cond( :,:,1 ) + log( model.a0 ), h.sector_log_pjoint_cond( :,:,2 ) + log( 1 - model.a0 ), h.sector_log_pjoint_cond( :,:,3 ) + log( 1 - model.a0 ) );
    a_mat = exp( [ permute( h.sector_log_pjoint, [1 3 2] ) ] );
    a_mat = [ a_mat  sum( a_mat, 2 ) ];
    a_mat = a_mat / sum( a_mat(2:end,end) ) * sum( h.w(2:end) );
    
    a_mat( 1,1 ) = a_mat( 1,1 ) + model.a0 * model.w00d;
    a_mat( 1,2 ) = a_mat( 1,2 ) + ( 1-model.a0 ) * model.w01d;
    
    figure; set( bar( h.x, h.w ), 'edgecolor', [ 0.7 0.7 0.7 ], 'facecolor', [ 0.7 0.7 0.7 ] );
    hold on; plot( h.x, a_mat( :,1:end-1), 'linewidth', 2 );
    plot( h.x, a_mat( :,end ), 'c--', 'linewidth', 1 );
    
    xlabel( 'ws.block.reduced\_ssm\_activeness(:)' );
    title( 'ws' )
    
    % Show the distributions in log domain
    
    [hlog.w, hlog.x] = hist( log( ws.block.reduced_ssm_activeness( find( ws.block.reduced_ssm_activeness > 0 ) ) ), ws.PAR.DETECTION.DEBUG_HIST_NBINS );
    hlog.x = hlog.x(:);
    hlog.w = running_mean( hlog.w, 3 );
    hlog.w = hlog.w(:) / sum( hlog.w );
    
    [ hlog.frame_log_post, hlog.sector_log_post, hlog.alldata_log_plkld, hlog.frame_log_plkld, hlog.frame_log_pjoint, hlog.sector_log_pjoint_cond ] = compute_post_2d2g1se( model, exp( hlog.x ) );
    hlog.sector_log_pjoint = cat( 3, hlog.sector_log_pjoint_cond( :,:,1 ) + log( model.a0 ), hlog.sector_log_pjoint_cond( :,:,2 ) + log( 1 - model.a0 ), hlog.sector_log_pjoint_cond( :,:,3 ) + log( 1 - model.a0 ) );
    a_mat = exp( repmat( hlog.x, 1, 3 ) + [ permute( hlog.sector_log_pjoint, [1 3 2] ) ] );
    a_mat = [ a_mat  sum( a_mat, 2 ) ];
    a_mat = a_mat / sum( a_mat(2:end,end) ) * sum( hlog.w(2:end) );
    
    a_mat( 1,1 ) = a_mat( 1,1 ) + model.a0 * model.w00d;
    a_mat( 1,2 ) = a_mat( 1,2 ) + ( 1-model.a0 ) * model.w01d;
    
    figure; set( bar( hlog.x, hlog.w ), 'edgecolor', [ 0.7 0.7 0.7 ], 'facecolor', [ 0.7 0.7 0.7 ] );
    hold on; plot( hlog.x, a_mat( :,1:end-1), 'linewidth', 2 );
    plot( hlog.x, a_mat( :,end ), 'c--', 'linewidth', 1 );
    
    xlabel( 'log( ws.block.reduced\_ssm\_activeness(:) )' );
    title( 'debug' )
    
    % Pause
    
    keyboard
    
  end
  
  
  
%--------------------------------------------------------------------------------
% Function used by EM in "fit_model" (above)
%
% Evaluate the negative of the tied term in log lkld of the whole
% data, that involves a Gamma + Shifted Rice mixture:
% - strictly positive data.
% - parameters gamma, beta, sigma and V.
%
% if isempty( ws.negloglike_w1 ), Gamma is fitted alone, without shifted Rice.

function y = negloglike( X )
  
  global ws;
  
  % Get the parameters
  
  t = X( 1 );

  gamma = X( 2 );
  beta  = t / gamma;
  
  beta_1  = X( 3 );
  
  % Impose constraint on them
  
  if ( gamma <= 1 ) | ( beta <= 0 ) | ( beta_1 <= 0 )
    y = +Inf; return;
  end
  
  % Do it!
  
  grt_ind = find( ws.block.reduced_ssm_activeness > t );
  
  ws.negloglike_w0 = ws.negloglike_w0(:);
  ws.negloglike_w1 = ws.negloglike_w1(:);
  
  y = 0;
  
  
  % Gamma
  if ~isempty( ws.negloglike_w0 )
    
    y = y - ws.negloglike_w0( ws.nz_ind(:) ).' * ( ( gamma - 1 ) .* ws.log_data_nz_ind - ws.block.reduced_ssm_activeness ( ws.nz_ind(:) ) / beta - gammaln( gamma ) - gamma * log( beta ) );
    
  end

  % Shifted Erlang
  
  if ~isempty( ws.negloglike_w1 )
    
    tmp = ws.block.reduced_ssm_activeness( grt_ind ) - t;
    
    y = y - ws.negloglike_w1( grt_ind(:) ).' * ( log( tmp(:) ) - tmp(:) / beta_1 - 2 * log( beta_1 ) );
    
  end
  
  % Normalize y by its number of terms,
  % so that "TolFun" may have a meaning
  
  y = y / numel( ws.block.reduced_ssm_activeness(:) );
  


  
% ----------------------------------------------------------------------
% Unsupervised detection threshold selection
% - one threshold for FAR = PAR.DETECTION.FAR_T
% - one threshold for FRR = PAR.DETECTION.FRR_T
% (IDIAP RR 05-52)

function select_thresholds
  
  global ws
  
  % Ordered posteriors (reduced data)
  
  [ log_post_act, sort_ind ] = sort( ws.block.reduced_ssm_log_post_act(:) );
  post_act = exp( log_post_act );
  post_sil = 1 - post_act;
  
  % Trying each sample as a strict threshold ( x > thr ),
  % we'll estimate the number of TPs, FPs, TNs, FNs
  
  N_TP = flipud( cumsum( flipud( [ post_act( 2:end ); 0 ] ) ) );
  N_FP = flipud( cumsum( flipud( [ post_sil( 2:end ); 0 ] ) ) );
  
  N_TN = cumsum( post_sil );
  N_FN = cumsum( post_act );
  
  % Derive the estimated FAR and estimated FRR
  
  [ far, frr ] = far_frr( N_TP, N_FP, N_TN, N_FN );
  
  % Pick the ones closest to our targets
  
  % FAR_T (we remove as much noise as possible)
  
  [ tmp, far_min_ind ] = min( abs( far - ws.PAR.DETECTION.FAR_T ) );
  
  ws.block.far_threshold = ws.block.reduced_ssm_log_post_act( sort_ind( far_min_ind ) );

  % FRR_T (we keep as much speech as possible)
  
  [ tmp, frr_min_ind ] = min( abs( frr - ws.PAR.DETECTION.FRR_T ) );
  
  ws.block.frr_threshold = ws.block.reduced_ssm_log_post_act( sort_ind( frr_min_ind ) );

  
  if ws.PAR.DETECTION.DEBUG
    
    figure; hist( exp( ws.block.reduced_ssm_log_post_act(:) ),100)
    
    figure; plot( far, frr, '.-' );
    hold on; plot( far( far_min_ind ), frr( far_min_ind ), 'r+', 'markersize', 15 );
    hold on; plot( far( frr_min_ind ), frr( frr_min_ind ), 'r+', 'markersize', 15 );
    
    
    disp( sprintf( 'selection: estimated FAR:%g%%, est. FRR:%g%%', 100 * far( far_min_ind ), 100 * frr( far_min_ind ) ) );
    disp( sprintf( 'selection: estimated FAR:%g%%, est. FRR:%g%%', 100 * far( frr_min_ind ), 100 * frr( frr_min_ind ) ) );
    
    keyboard
    
  end
  

% -----------------------------------------------------------------------
% Extract azimuth differences from the last block,
% that will be used to fit a bi-Gaussian model for short-term clustering.
% (see "reduce_stc_data" and "fit_stc_model")

function extract_stc_data
  
  global ws
  
  % ( 1 ) Extract (time, azimuth) values.

  t  = ws.block.sphloc_mat( 1, : );
  az = ws.block.sphloc_mat( 3, : );
  
  % ( 2 ) Compare them, for each possible time difference
  
  dt_mat  = abs( repmat(  t(:), 1, numel( t ) ) - repmat(  t(:).', numel( t ), 1 ) );
  % Don't compare a location estimate with itself
  % Don't compare twice the same pair of location estimates
  dt_mat( find( tril( ones( numel( t ) ) ) ) ) = NaN;
  
  daz_mat = repmat( az(:), 1, numel( t ) ) - repmat( az(:).', numel( t ), 1 );

  delta_t_max = min( max( abs( dt_mat(:) ) ), ws.PAR.STC.T_SHORT );
  ws.block.stc_data = cell( delta_t_max, 1 );
  
  for delta_t = 0:delta_t_max

    ind = find( dt_mat == delta_t );
    
    ws.block.stc_data{ delta_t + 1 } = daz_mat( ind );

  end
  


% ----------------------------------------------------------------------
% Fading memory applied to the azimuth differences
% given by "extract_stc_data"

function reduce_stc_data

  global ws

  n_previous = numel( ws.block.previous_stc_data );
  n_current  = numel( ws.block.stc_data );
  n_max      = max( n_previous, n_current );

  % ( 0 ) initialize

  ws.block.reduced_stc_data = cell( n_max, 1 );

  for a = 1:n_max
    
    ws.block.reduced_stc_data{ a } = zeros( 1, 0 );

  end

  % Note: the code is replicated to avoid data replication
    
  % ( 1 ) Reduce "previous data" and copy it
  
  for a = 1:n_previous

    n_loc_est     = size( ws.block.previous_stc_data{ a }, 2 );
    max_n_loc_est = ws.a_struct.datarepr_maxprevious_frames; 
    
    if n_loc_est <= max_n_loc_est
    
      % Nothing to do
      
      ws.block.reduced_stc_data{ a } = ws.block.previous_stc_data{ a };
      
    else
      
      % Reduction is necessary: drop some frames, at regular time intervals
      
      ind = min( n_loc_est, 1 + round( ( n_loc_est-1 ) * ( 0.5:max_n_loc_est ) / max_n_loc_est ) );
      
      ws.block.reduced_stc_data{ a } = ws.block.previous_stc_data{ a }( :, ind );
      
    end
  
  end
  
  % ( 2 ) Reduce "current data" and append it
  
  for a = 1:n_current
    
    n_loc_est     = size( ws.block.stc_data{ a }, 2 );
    max_n_loc_est = ws.a_struct.datarepr_maxcurrent_frames; 
    
    if n_loc_est <= max_n_loc_est
      
      % Nothing to do
      
      ws.block.reduced_stc_data{ a } = [ ws.block.reduced_stc_data{ a }(:);   ws.block.stc_data{ a }(:)  ];
    
    else
    
      % Reduction is necessary
      
      ind = min( n_loc_est, 1 + round( ( n_loc_est-1 ) * ( 0.5:max_n_loc_est ) / max_n_loc_est ) );
      
      tmp = ws.block.stc_data{ a }( ind );

      ws.block.reduced_stc_data{ a } = [ ws.block.reduced_stc_data{ a }(:);  tmp(:) ];
      
    end

  end
  
  % Sanity check

  for a = 1:n_max

    if size( ws.block.reduced_stc_data{ a }, 2 ) > ws.a_struct.datarepr_maxprevious_frames + ws.a_struct.datarepr_maxcurrent_frames

      error( [ mfilename '.reduce_stc_data(): insanity detected.' ] );

    end

  end

    


% ----------------------------------------------------------------------
% Conceptually simple, but annoying to implement:
% bi-Gaussian EM training + constraints described in IDIAP RR 04-14
%
% This function modifies ws.stc.model
%
% Special note for delta_t = 0: since multiple simultaneous estimates of the same
% source are rare, we don't train a model for delta_t = 0, but simply re-use the
% model trained for delta_t = 1.

function update_stc_model

  global ws
  
  % We start from scratch every time, only the fading memory implemented
  % in "reduce_stc_data" provides us with some constraint from the past.
  
  % Update only if we have data
  
  nz_delta_t = [];
  
  delta_t_max = min( numel( ws.block.reduced_stc_data ) - 1, ws.PAR.STC.T_SHORT );
  
  % Ignore delta_t = 0 for now.
  
  for delta_t = 1:delta_t_max
    
    if numel( ws.block.reduced_stc_data{ delta_t + 1 } ) > 2

      nz_delta_t( end+1 ) = delta_t;

    end

  end

  % ( 1 ) Initialization of the EM algorithm
  % The main component is azimuth, used to try various thresholds.
  % Only azimuth is used to evaluate lkld ( as in step 2 ).
  
  n1 = max( 1, round( ws.PAR.STC.AUTO_INIT_N / 2 ) );
  n2 = max( 1, ws.PAR.STC.AUTO_INIT_N - n1 );

  % For each delta_t
  
  for delta_t = nz_delta_t(:).'
    
    some_data = ws.block.reduced_stc_data{ delta_t + 1 };
    some_data = some_data(:).';

    % Sort the data by increasing abs( delta azimuth )
    
    [ v, ind ] = sort( abs( some_data ) );
    
    some_data = some_data( ind );
    
    % Prepare a list of tentative thresholds on azimuth
    
    v = v( find( v < ws.PAR.STC.DAZ_MAX ) );
    
    % Half of them at equal intervals between min and max
    
    tmp_n1 = min( numel( v ) - 1, n1 );
    
    thr_list = v( 1 ) + diff( v( [ 1 end ] ) ) * ( 1:tmp_n1 ) / ( tmp_n1 + 1 );
    
    % The other half at equal percentile intervals 
    
    tmp_n2 = min( numel( v ) - 1, n2 );
    
    n = numel( v );
    ind = min( n, max( 1, 1 + ( n - 1 ) * ( 0.5:tmp_n2 ) / tmp_n2 ) );
    thr_list = unique( [ thr_list  interp1( 1:n, v, ind, 'linear' ) ] );
    
    % Try each of them

    n_thr = numel( thr_list );
    
    thr_result = [];
    
    for i_thr = 1:n_thr
      
      thr_result( i_thr ).thr = thr_list( i_thr );
      
      is_same = ( v < thr_result( i_thr ).thr );
      
      % priors
      
      thr_result( i_thr ).p_same    = mean( is_same );
      thr_result( i_thr ).p_notsure = 1 - thr_result( i_thr ).p_same;
      
      thr_result( i_thr ).alldata_log_lkld = -Inf;
      
      if ( thr_result( i_thr ).p_same > 0 ) & ( thr_result( i_thr ).p_notsure > 0 )
	
	some_data_same    = some_data( :,  is_same );
	some_data_notsure = some_data( :, ~is_same );
	
	% parameters
	  
	thr_result( i_thr ).same    = max( ws.PAR.STC.MODEL_STD_FLOOR, sqrt( mean( some_data_same .^ 2,    2 ) ) );
	thr_result( i_thr ).notsure = max( ws.PAR.STC.MODEL_STD_FLOOR, sqrt( mean( some_data_notsure .^ 2, 2 ) ) );
	  
	% Verify that the constraint imposed in IDIAP RR 04-14 is verified
	
	if all( thr_result( i_thr ).notsure  >  thr_result( i_thr ).same )
	  
	  % lkld of the whole data
	  
	  old_warning_state = warning;
	  warning off;
	  
	  tmp1 = log( thr_result( i_thr ).p_same )    + log_normpdf( some_data.', 0, thr_result( i_thr ).same );
	  tmp2 = log( thr_result( i_thr ).p_notsure ) + log_normpdf( some_data.', 0, thr_result( i_thr ).notsure );
	  
	  warning( old_warning_state );
	  
	  thr_result( i_thr ).alldata_log_lkld = sum( my_logsum_fast( tmp1, tmp2 ) );

	end

      end
      
    end  % of loop through thresholds
      
      
    [ tmp, best_ind ] = max( [ thr_result.alldata_log_lkld ] );
    
    if tmp > -Inf
      
      % Usual case: initialization worked for at least one threshold
      % -> pick the best one in terms of likelihood and update the model.
      
      ws.stc.model.same( :, delta_t+1 )       = thr_result( best_ind ).same;
      ws.stc.model.notsure( :, delta_t+1 )    = thr_result( best_ind ).notsure;
      
      ws.stc.model.p_same( delta_t+1 )     = thr_result( best_ind ).p_same;
      ws.stc.model.p_notsure( delta_t+1 )  = thr_result( best_ind ).p_notsure;

    end
    
    
  end  % of loop through "delta_t"
  


  
  % ( 2 ) EM algorithm

  % For each delta_t

  for delta_t = nz_delta_t(:).'
    
    % Get the azimuth differences
    
    x = ws.block.reduced_stc_data{ delta_t + 1 };
    x = x(:).';

    old_warning_state = warning;
    warning off;  
    two_log_abs_x = 2 * log( abs( x ) ).';
    warning( old_warning_state );
    
    % Get the initialized parameters
    
    p_same    = ws.stc.model.p_same( delta_t + 1 );
    p_notsure = ws.stc.model.p_notsure( delta_t + 1 );
    
    same    = ws.stc.model.same( :, delta_t + 1 );
    notsure = ws.stc.model.notsure( :, delta_t + 1 );
    
    % where we'll store the new parameters
    % when we are sure that lkld increases
    %
    % Indeed, since we are forcing means to zero,
    % global lkld is not guaranteed to increase.

    new_p_same = p_same;
    new_p_notsure = p_notsure;
    
    new_same   = same;
    new_notsure = notsure;

    % For the "ws.PAR.DAZ_MAX" condition
    
    ind_out_of_bound = find( abs( x ) > ws.PAR.STC.DAZ_MAX );
    
    % Start EM
    
    previous_alldata_log_lkld = -Inf;
    iter = 0;
    finished = 0;
    while ~finished
      
      % ( 2.1 ) E step: compute posteriors

      old_warning_state = warning;
      warning off;
      log_joint_same    = log( p_same )    + log_normpdf( x.', 0, same );
      log_joint_notsure = log( p_notsure ) + log_normpdf( x.', 0, notsure );
      warning( old_warning_state );
      
      % Impose the "ws.PAR.DAZ_MAX" condition in a hard manner
      log_joint_same( ind_out_of_bound ) = -Inf;
      
      log_lkld = my_logsum_fast( log_joint_same, log_joint_notsure );
      alldata_log_lkld = sum( log_lkld );

      log_post_same    = log_joint_same    - log_lkld;
      log_post_notsure = log_joint_notsure - log_lkld;

      % Sanity check
      
      tmp = my_logsum_fast( log_post_same, log_post_notsure );
      if any( abs( tmp ) > 1e-8 )
	error( 'post must sum to 1' );
      end

      % If we have an increase, let's update the parameters

      if alldata_log_lkld > previous_alldata_log_lkld

	new_p_same    = p_same;
	new_p_notsure = p_notsure;
	new_same      = same;
	new_notsure   = notsure;

      end
      

      % ( 2.2 ) M step: update parameters
      
      log_total = [ my_logsum_fast( log_post_same )   my_logsum_fast( log_post_notsure ) ];
      tmp = log_total - my_logsum_fast( log_total );
      tmp = exp( tmp );
      
      % Correct minor numerical errors
      tmp = tmp / sum( tmp );
      
      % Update priors
      p_same    = tmp( 1 );
      p_notsure = tmp( 2 );
      
      % Update standard deviations

      tmp = log_post_same(:) + two_log_abs_x;
      same = exp( 0.5 * (  my_logsum_fast( tmp ) -  log_total( 1 )  ) );

      tmp = log_post_notsure(:) + two_log_abs_x;
      notsure = exp( 0.5 * (  my_logsum_fast( tmp ) -  log_total( 2 )  ) );
      
      % ( 2.3 ) Termination test
    
      iter = iter + 1;
      finished = ( iter >= ws.PAR.STC.UPDATE_LOCAL_MAX_ITER ) | ( abs( alldata_log_lkld - previous_alldata_log_lkld ) / abs( alldata_log_lkld ) < ws.PAR.STC.UPDATE_LOCAL_CV_THR );
      
      previous_alldata_log_lkld = alldata_log_lkld;
      
    end  % of EM iterations
    
    % Set the parameters

    ws.stc.model.p_same( delta_t + 1 )     = new_p_same;
    ws.stc.model.p_notsure( delta_t + 1 )  = new_p_notsure;

    ws.stc.model.same( :, delta_t + 1 )    = new_same;
    ws.stc.model.notsure( :, delta_t + 1 ) = new_notsure;

  end    % of loop through delta_t's


  % ( 3 ) Enforce   same < notsure  condition
  % ( not supposed to be necessary, given that we already check it at initialization )

  ws.stc.model.notsure = max( ws.stc.model.notsure, 1.01 * ws.stc.model.same );

  % ( 4 ) Model for delta_t = 0
  
  ws.stc.model.p_same( 1 )    = mean( ws.stc.model.p_same( 2:end ) );
  ws.stc.model.p_notsure( 1 ) = 1 - ws.stc.model.p_same( 1 );

  ws.stc.model.same( 1 )    = mean( ws.stc.model.same( 2:end ) );
  ws.stc.model.notsure( 1 ) = mean( ws.stc.model.notsure( 2:end ) );

  

% ----------------------------------------------------------------------
% Apply binary decision (SAM-SPARSE-MEAN) in the spectral domain.
% -> for each sector, a different binary mask.
% 
% Input:
%   ws.block.X
%   ws.block.ssm_binary_decision
% 
% Output:
%   ws.block.abs_X_sector


function extract_abs_X_sector

  global ws

  ws.block.abs_X_sector = repmat( NaN, [ ws.PAR.FFT_FRAME_LENGTH, ws.a_struct.n_sectors, ws.block.n_frames ] );

  abs_X = abs( ws.block.X( :, ws.PAR.MFCC_CHANNEL, : ) );

  for a = 1:ws.a_struct.n_sectors
    
    ws.block.abs_X_sector( :, a, : ) = abs_X .* double( ws.block.ssm_binary_decision( :, a, : ) );
    
  end
  
  % Sanity check
  if any( isnan( ws.block.abs_X_sector ) )
    error( [ mfilename '.extract_abs_X_sector() is NaN-insane.' ] );
  end

  % Sanity check
  if any( isinf( ws.block.abs_X_sector ) )
    error( [ mfilename '.extract_abs_X_sector() is Inf-insane.' ] );
  end

  


% ----------------------------------------------------------------------
% Update data representation for USS
%
% Input: 
%   ws.block.X
%
% Output:
%   ws.block.datarepr

function update_datarepr_uss

  global ws

  % "1uss" -> only one model, using the original spectrum abs( ws.block.X( :, ws.PAR.MFCC_CHANNEL, : ) )
  % of this microphone. This way we have a chance to capture the background noise.
  
  % Reduce past data
  
  old_data = [];
  if isfield( ws.block, 'datarepr' )
    old_data = ws.block.data_repr;
  end
  
  n = numel( old_data );
  n_max = ws.PAR.USS.data_par.past_npoints_max;
  
  if n > n_max
    
    old_data = sort( old_data );
    
    old_data = old_data( 1 + round( ( n-1 ) * ( 0:n_max-1 ) / ( n_max - 1 ) ) );
    
    % The distribution stays approximately the same
    
  end
  
  % Reduce new data
  
  new_data = abs( ws.block.X( :, ws.PAR.MFCC_CHANNEL, : ) );
  new_data = new_data( find( new_data ) );      % non-zero only
  
  n = numel( new_data );
  n_max = ws.PAR.USS.data_par.future_npoints;
  
  if n > n_max
    
    new_data = sort( new_data );
    
    new_data = new_data(  1 + round( ( n-1 ) * ( 0:n_max-1 ) / ( n_max - 1 ) ) );
    
    % The distribution stays approximately the same
    
  end
  
  % Merge the two and update
  
  ws.block.data_repr = [ old_data(:); new_data(:) ];
  
  
% ----------------------------------------------------------------------
% Unsupervised Spectral Subtraction (USS)
%
% Inputs:
%   ws.block.abs_X_sector
%   ws.block.datarepr
%
% Output:
%   ws.block.raylsherl
%   ws.block.abs_X_sector (modified)

function remove_noise

  global ws

  % A single model on the original spectrum
  % ( this way we can observe the background noise )
  
  if numel( ws.block.data_repr ) > 3
    
    % Fit a Rayleigh + Shifted Erlang model
    
    ws.block.raylsherl = fit_raylsherl( ws.block.data_repr, ws.PAR.USS.em_par );
    
  else
    
    if ws.PAR.USS.em_par.verbose
      disp( sprintf( [ mfilename '.remove_noise(): zero data in a block, on frames %d-%d.' ], ws.block.start_frame, ws.block.stop_frame ) );
    end
    
    % Not enough non-zero samples to fit a model
    
    % Most of the data is zero => we consider this as noise
    
    tmp.erlang_t = +Inf;
    ws.block.raylsherl = tmp;
    
  end
  
  % Now we can apply this unique model
  % to all sector spectra

  for a = 1:ws.a_struct.n_sectors
    
    % Spectral "subtraction"
    % ( not exactly subtraction, but the best choice we found
    % for subsequent cepstrum extraction ).
    
    ws.block.abs_X_sector( :,a,: ) = max( 1, ws.block.abs_X_sector( :,a,: ) / ws.block.raylsherl.erlang_t );
    
  end
  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix to linearly transform magnitude values
% into MEL filterbank values

function fbank_mat = prepare_fbank_mat( par )

  if nargin < 1
    error( [ mfilename ' needs 1 input argument.' ] );
  end

  % Check for mandatory parameters
  required_fields = { 'fs', 'fft_framelength', 'mfcc_band_hz', 'mfcc_B' };
  check_param( required_fields, fieldnames( par ) );
  
  % Define the MEL scale, following HTK book 3.3

  mel_from_f = inline( '2595 * log10( 1 + f / 700 )', 'f' );
  f_from_mel = inline( '700 * ( 10 .^ ( mel / 2595 ) - 1 )', 'mel' );

  % Sanity check
  f = 10000;
  if abs( f_from_mel( mel_from_f( f ) ) - f ) > 1e-8
    error( [ mfilename '.prepare_fbank_mat() is insane.' ] );
  end
  
  melband = mel_from_f( max( 0, min( par.fs / 2, par.mfcc_band_hz ) ) );
  
  % Equal spacing in MEL space along the band

  n = par.mfcc_B + 1;
  mel_list = melband( 1 ) + diff( melband ) * ( 0:n ) / n;
  
  % Back to frequency scale

  f_list = max( par.mfcc_band_hz( 1 ), min( par.mfcc_band_hz( 2 ), f_from_mel( mel_list ) ) );
  
  % Define the triangular filters
  
  fbank_mat = zeros( par.mfcc_B, 2 * par.fft_framelength );

  f_k = ( 0:( 2 * par.fft_framelength-1 ) ) / ( 2 * par.fft_framelength ) * par.fs;
  
  for a = 1:par.mfcc_B
    
    [ tmp, k_low ]    = min( abs( f_k - f_list( a ) ) );
    [ tmp, k_center ] = min( abs( f_k - f_list( a + 1 ) ) );
    [ tmp, k_high ]   = min( abs( f_k - f_list( a + 2 ) ) );

    fbank_mat( a, k_low:k_high ) = interp1( [ k_low  k_center  k_high ], [ 0 1 0 ] / max( 1, k_high - k_low ), k_low:k_high, 'linear', NaN );

  end
  
  % Sanity check
  if any( isnan( fbank_mat(:) ) )
    error( [ mfilename '.prepare_fbank_mat() is still insane.' ] );
  end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrix to linearly transform log filterbank values
% into cepstral coefficients

function dct_mat   = prepare_dct_mat( par )

  if nargin < 1
    error( [ mfilename '.prepare_dct_mat() needs 1 input argument.' ] );
  end

  % Check for mandatory parameters
  required_fields = { 'mfcc_B', 'mfcc_N', 'mfcc_L' };
  check_param( required_fields, fieldnames( par ) );
  
  % DCT: following formula in HTK book 3.3
  
  r = 0:par.mfcc_N;  % Includes the C_0 coefficient
  c = 1:par.mfcc_B;
  
  dct_mat = sqrt( 2 / par.mfcc_B ) * cos( pi / par.mfcc_B * r(:) * ( c(:).' - 0.5 ) );
  
  % Optional liftering: following formula in HTK book 3.3
  
  if ~isempty( par.mfcc_L )

    if par.mfcc_L ~= 0
      
      dct_mat = full( diag( sparse( 1 + par.mfcc_L / 2 * sin( pi * r / par.mfcc_L ) ) ) * dct_mat );
      
    end
    
  end
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Implementation of MFCC computation

function [ mfcc, fbank ] = compute_mfcc( mag_spectrum )

  if nargin < 1
    error( [ mfilename '.compute_mfcc() needs 1 input argument.' ] );
  end

  if isempty( mag_spectrum )

    fbank = zeros( size( ws.a_struct.fbank_mat, 1 ), 0 );
    mfcc  = zeros( size( ws.a_struct.dct_mat, 1 ), 0 );
    return;
    
  end

  % "ws" only used to access parameters (in ws.PAR),
  % as well as precomputed transformation matrices (in ws.a_struct)
  % (filterbanks and DCT).
  %
  % All the stuff in "ws.block" is NOT modified at all.

  global ws

  % From linear spectrogram to filterbank spectrogram
  
  if ~ws.PAR.MFCC.usepower
    
    fbank = ws.a_struct.fbank_mat * mag_spectrum;
    
  else

    fbank = ws.a_struct.fbank_mat * ( mag_spectrum .^ 2 );
    
  end

  % Update minimum observed non-zero value
  
  tmp = fbank;
  tmp( find( tmp <= 0 ) ) = +Inf;
  tmp = min( tmp, [], 2 );
  
  tmp( find( isnan( tmp ) | isinf( tmp ) | ( tmp <= 0 ) ) ) = ws.PAR.MFCC.fbank_floor;
  
  if ~isfield( ws, 'previous_min_fbank_mag' )

    ws.previous_min_fbank_mag = [];

  end

  if isempty( ws.previous_min_fbank_mag )
    
    ws.previous_min_fbank_mag = tmp;
    
  else
    
    ws.previous_min_fbank_mag = min( ws.previous_min_fbank_mag, tmp );
    
  end
  
  % Sanity check
  if any( ws.previous_min_fbank_mag <= 0 )
    error( [ mfilename '.compute_mfcc() is insane!' ] );
  end

  % Compute MFCCs
  
  mfcc = ws.a_struct.dct_mat * log( max( repmat( ws.previous_min_fbank_mag, 1, size( fbank, 2 ) ), fbank ) );
  