function ssm_result = extract_ssm_fast( wavefilename_template, ssm_parameters_filename, block_size_frames, timeframe_overlap )
  
% function ssm_result = extract_ssm_fast( wavefilename_template, ssm_parameters_filename,  block_size_frames, timeframe_overlap )
%
% block_size_frames is an optional parameter: number of frames to process in a block (default 1000)
%
% 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( 'block_size_frames', 'var' )
    block_size_frames = 1000;
  end
  
  if ~exist( 'timeframe_overlap', 'var' )
    timeframe_overlap = 0.5;
  end
  

  % To pass large pieces of data efficiently,
  % to and from the various subfonctions
  % ( namely the C implementations of GCC-PHAT and SAM-SPARSE-MEAN )
  global ws  % Stands for "workspace"

  ws = [];

  
  % ( 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.block_size_frames       = block_size_frames;
  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 );

  % Prepare parameters for the GCC-PHAT C implementation
  
  ws.a_struct.wavefile_list = ssm_result.wavefile_fullname_list;
  ws.a_struct.n_channels    = numel( ssm_result.wavefile_fullname_list );
  ws.a_struct.framelength   = ssm_result.framelength;
  ws.a_struct.frameshift    = ssm_result.frameshift;

  ws.a_struct.win           = diag( sparse( hamming( ws.a_struct.framelength, 'periodic' ) ) );
  
  ws.a_struct.block_size_frames = block_size_frames;  
  ws.a_struct.start_frame   = 1;
  ws.a_struct.stop_frame    = nframes;
  
  ws.PAR.GCCPHAT_DO_ZERO_MEAN = 1;
  ws.PAR.GCCPHAT_PREEMP     = PREEMP;
  
  ws.a_struct.n_pairs       = size( ssm.pairs, 1 );
  
  % Prepare parameters for the SAM-SPARSE-MEAN C implementation
  
  ws.a_struct.n_sectors     = numel( ssm.sector );
  ws.a_struct.conj_ssm_Z    = conj_ssm_Z;
  ws.ssm = ssm;
  ws.ssm.pairs = uint16( ws.ssm.pairs );
  
  % Count time
  
  
  ssm_result.total_cputime = 0;
  ssm_result.total_gccphat_cputime = 0;
  ssm_result.total_ssm_cputime = 0;
  
  start_cputime = cputime;
    
  % Start iterating through frames, processing them in a blockwise fashion
  
  if VERBOSE
    chrono = chrono_start( nframes, 30 );
  end
  
  for i_frame = ws.a_struct.start_frame:ws.a_struct.stop_frame
    
    % We do as much processing as possible per block of frames.
    is_new_block = ~mod( i_frame - 1, block_size_frames );
    
    % Absolute frame index from the beginning of the file (starting with 1)
    frame_ind    = i_frame;
    

    % ( 4.1 ) Get the complex-domain GCC-PHAT measurements
    % ( C implementation )

    if is_new_block
      
      ws.block.start_frame = frame_ind;
      ws.block.stop_frame  = min( ws.block.start_frame + ws.a_struct.block_size_frames - 1, ws.a_struct.stop_frame );
      
      if VERBOSE
	disp( ' ' );
	disp( sprintf( [ 'Extracting GCC-PHAT from waveforms (frames %d to %d)...' ], ws.block.start_frame, ws.block.stop_frame ) );
      end
      
      get_observed_gccphat;
 
    end

    
        
    % ( 4.2 ) Extract the SSM activeness
    % ( C implementation )
    
    if is_new_block
      
      if VERBOSE
	disp( sprintf( [ 'Extracting SSM activeness from GCC-PHAT (frames %d to %d)...' ], ws.block.start_frame, ws.block.stop_frame ) );
      end
      
      get_ssm_activeness;
            
    end

    % ( 4.3 ) Store the result

    ssm_result.activeness( :, ws.block.start_frame:ws.block.stop_frame ) = ws.block.ssm_activeness;
    
    if VERBOSE
      chrono = chrono_check( chrono, i_frame );
    end
    
  end    
  
  if VERBOSE
    chrono_stop( chrono );
    disp( ' ' );
    disp( [ mfilename ' is done.' ] );
    disp( ' ' );
  end
  
  
  
% ----------------------------------------------------------------------
% Wrapper for the C implementation of GCC-PHAT
% Get the complex-domain GCC-PHAT measurements

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;
  
  % 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.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.GCCPHAT_DO_ZERO_MEAN
      sig = sig - repmat( mean( sig, 1 ), size( sig, 1 ), 1 );
    end
    
    % Pre-emphasis
    sig = sig( 2:end, : ) - ws.PAR.GCCPHAT_PREEMP * sig( 1:end-1, : );
    
    % zero-padded FFT of hamming-windowed signal
    tmp = fft( ws.a_struct.win * sig, 2 * ws.a_struct.framelength, 1 );
    ws.block.X( :,i_channel,: ) = permute( tmp( 2:ws.a_struct.framelength+1, : ), [ 1 3 2 ] );
    
  end
  
  clear sig tmp

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


  if 1

    % C implementation (26 times faster and less memory-consuming)
    s051115_compute_gccphat;

  else

    % Original MATLAB implementation
    % Not the smartest way, but I'm trying to minimize memory usage
    ws.block.ssm_observed_gccphat = exp( j * angle( ws.block.X( :, ws.ssm.pairs( :,1 ), : ) .* conj( ws.block.X( :, ws.ssm.pairs( :,2 ), : ) ) ) );

  end
  
  % Clear out memory
  ws.block.X = [];
  
  
  
  
    
% ----------------------------------------------------------------------
% Wrapper for the C implementation of SAM-SPARSE-MEAN
%
% 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" ).

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.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.

    s051114_get_ssm_activeness;

  else


    % Original MATLAB implementation

    % To do it in one step, we would need a large amount of memory
    % ( 4-D )
    % -> loop through frames
    % -> eventually we'll do a C implementation
    
    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
