function ssm_result = extract_ssm( wavefilename_template, ssm_parameters_filename, timeframe_overlap )
  
% function ssm_result = extract_ssm( wavefilename_template, ssm_parameters_filename, timeframe_overlap )
%
% timeframe_overlap is an optional parameter between 0 and 1 (default value 0.5, which means 50% overlap
% between consecutive time frames).
%
% Return a ssm_result structure with the following fields:
% ssm_result.wavefilename_template       string
% ssm_result.wavefile_fullname_list      cell array of NCHANNELS string 
% ssm_result.ssm_parameters_filename     string
% ssm_result.timeframe_overlap           scalar (0 to 1)
% ssm_result.nsamples                    integer
% ssm_result.fs                          integer
% ssm_result.framelength                 integer
% ssm_result.frameshift                  integer
% ssm_result.total_duration              scalar (seconds)
% ssm_result.activeness                  NSECTORS x NFRAMES matrix of activeness values (between 0 and 1).
%
% This script takes the wave files matching "wavefilename_template", orders them alphanumerically,
% loads the SAM-SPARSE-MEAN parameters in "ssm_parameters_filename".
%
% SAM-SPARSE-MEAN "activeness" values are then extracted for each time frame and 
% each sector of space, and stored into the "ssm_result.activeness" matrix.
%
% Each value is between 0 and 1, indicating in how much of the bandwidth (e.g. 0.75 * 8 kHz )
% acoustic sources in a given sector of space are dominant.
%
% For further details on SAM-SPARSE-MEAN, see IDIAP RR 04-67.
%
% Guillaume Lathoud
% 2005

  if nargin < 2
    error( [ mfilename ' needs at least 2 input arguments.' ] );
  end
  
  dbstop error
  
  VERBOSE = 1;
  
  % Pre-emphasis DOES help, even if magnitude is not involved.
  PREEMP  = 0.97;
  
  if ~exist( 'timeframe_overlap', 'var' )
    timeframe_overlap = 0.5;
  end
  
  % ( 1 ) Find the wavefilenames
  
  d = dir( wavefilename_template );
  d = d( find( ~[ d.isdir ] ) );
  [ tmp, ind ] = sort( { d.name } );
  
  a_path = fileparts( wavefilename_template );
  

  if VERBOSE
    disp( ' ' );
    disp( [ mfilename ' found the following wave files (in this order):' ] );
  end
    
  
  wavefile_fullname_list = {};
  nsamples_list          = [];
  fs_list                = [];
  for a = 1:numel( d )
    
    % Build the complete filename
    
    wavefile_fullname_list{ a } = fullfile( a_path, d( a ).name );

    % Extract some information
    
    [ a_size, a_fs ] = wavread( wavefile_fullname_list{ a }, 'size' );
    
    nsamples_list( a ) = a_size( 1 );
    
    if a_size( 2 ) ~= 1
      error( [ mfilename ' requires single channel files only.' ] );
    end
    
    fs_list( a )       = a_fs;
    
    
    if VERBOSE
      disp( sprintf( '%d: "%s"', a, wavefile_fullname_list{ a } ) );
    end
      
  end
  
  % Check consistency of the sampling frequency
  
  fs = unique( fs_list );
  if numel( fs ) > 1
    error( [ mfilename ': the sampling frequency must be the same for all wave files.' ] );
  end
  
  % ( 2 ) Find and load the SSM parameters
  
  load( ssm_parameters_filename, 'ssm' );
  
  if VERBOSE
    disp( ' ' );
    disp( [ mfilename ' loaded the SSM parameter file "' ssm_parameters_filename '".' ] );
  end
  
  
  % ( 3 ) Consistency checks + find the total length
  
  if fs ~= ssm.fs
    error( [ mfilename ': the sampling frequency of wave files is not consistent with the SSM parameters in "' ssm_parameters_filename '".' ] );
  end
  
  channel_list = unique( ssm.pairs(:) );
  
  if ( min( channel_list ) < 1 ) | ( max( channel_list ) > numel( wavefile_fullname_list ) )
    error( [ mfilename ': inconsistent number of channels.' ] );
  end
  
  nchannels = max( ssm.pairs(:) );
  
  % Total length
  
  nsamples = min( nsamples_list );
  
  total_duration = nsamples / fs;
  
  if VERBOSE
    
    if numel( unique( nsamples_list ) ) > 1
      disp( sprintf( [ mfilename ': WARNING! Wave file do not have the same duration. Had to troncate down to the shortest duration.' ] ) );
    end
  
  end
  
  % Frame shift
  
  frameshift = round( ssm.framelength * ( 1 - timeframe_overlap ) );


  
  % Store parameters
  
  ssm_result.wavefilename_template   = wavefilename_template;
  ssm_result.wavefile_fullname_list  = wavefile_fullname_list;
  ssm_result.ssm_parameters_filename = ssm_parameters_filename;
  ssm_result.timeframe_overlap       = timeframe_overlap;
  ssm_result.nsamples                = nsamples;
  ssm_result.fs                      = fs;
  ssm_result.framelength             = ssm.framelength;
  ssm_result.frameshift              = frameshift;
  ssm_result.total_duration          = total_duration;
  ssm_result.win                     = diag( sparse( hamming( ssm_result.framelength, 'periodic' ) ) );
  
  if VERBOSE
    disp( '  ');
    disp( sprintf( [ mfilename ' parameters:' ] ) );
    disp( ssm_result );
  end
  
  if ssm_result.nsamples < ssm_result.framelength 
    error( [ mfilename ' does not have enough samples for even one frame.' ] );
  end
  
  
  % ( 4 ) Extract activeness through SAM-SPARSE-MEAN analysis

  conj_ssm_Z = conj( ssm.Z );
  
  nsectors = numel( ssm.sector );
  nframes  = 1 + floor( ( ssm_result.nsamples - ssm_result.framelength ) / ssm_result.frameshift );
  
  ssm_result.activeness = repmat( NaN, nsectors, nframes );

  
  
  % Start iterating through frames
  
  if VERBOSE
    chrono = chrono_start( nframes, 30 );
  end
  
  for i_frame = 1:nframes

    % Add one sample for pre-emphasis
    firstsample = max( 1, ( i_frame - 1 ) * ssm_result.frameshift );
    
    lastsample  = ssm_result.framelength + ( i_frame - 1 ) * ssm_result.frameshift;
    
    % ( 4.1 ) Load waveforms
    
    % "1+" for pre-emphasis
    
    x_mat = repmat( NaN, 1 + ssm_result.framelength, nchannels );
    
    for i_channel = channel_list(:).'
      
      x = wavread( ssm_result.wavefile_fullname_list{ i_channel }, [ firstsample  lastsample ] );
      
      if numel( x ) < ssm_result.framelength + 1
	% For the first frame: prepend a zero
	x = [ 0; x ];
      end
      
      x_mat( :,i_channel ) = x;
      
    end
    
    % Pre-emphasis
    
    x_mat = x_mat( 2:end, : ) - PREEMP * x_mat( 1:end-1, : );
    
    % Zero-mean
    
    x_mat = x_mat - repmat( mean( x_mat, 1 ), size( x_mat, 1 ), 1 );
    
    
    ssm_result.activeness( :, i_frame ) = 0;
    
    if sum( any( x_mat( :, channel_list ) ) ) > 1
    
      % ( 4.2 ) Compute GCC-PHAT
      
      % zero-padded, hamming-windowed FFT
      X_mat = fft( ssm_result.win * x_mat, 2 * ssm_result.framelength, 1 );
      
      % GCC-PHAT
      CC = X_mat( :, ssm.pairs( :,1 ) ) .* conj( X_mat( :, ssm.pairs( :,2 ) ) );

      GCCPHAT = zeros( size( CC ) );
      
      abs_CC  = abs( CC );
      nonzero = find( abs_CC > 0 );
      GCCPHAT( nonzero ) = CC( nonzero ) ./ abs_CC( nonzero );
      
      % ( 4.3 ) Compute SAM-SPARSE-MEAN activeness
      
      % We only need strictly positive frequencies
      
      GCCPHAT_posit = GCCPHAT( 2:ssm_result.framelength + 1, : );
      
      % "MEAN": Average delay-sum power in each sector,
      % at each frequency
      
      adsp = permute( mean( 0.5 * ( 1 + real( repmat( GCCPHAT_posit, [ 1 1 nsectors ] ) .* 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, nsectors ) );
      
      % If two sectors have same "adsp" value for a frequency bin, 
      % we ignore that frequency ( rare case ).
      
      ind_violation = find( sum( binary_decision, 2 ) > 1 );
      binary_decision( ind_violation, : ) = 0;
      
      % Activeness
      ssm_result.activeness( :, i_frame ) = mean( binary_decision, 1 ).';
      
    end

    if VERBOSE
      chrono = chrono_check( chrono, i_frame );
    end
    
  end    
  
  if VERBOSE
    chrono_stop( chrono );
    disp( ' ' );
    disp( [ mfilename ' is done.' ] );
    disp( ' ' );
  end