function out = powerdecay_analysis( p )

% function out = powerdecay_analysis( p )
% 
% Compute the fullband and narrowband power decay between lapel
% microphone and distant microphone. This is done at several
% locations.  Throughout this whole function, "location" means
% "spatial location".
% 
% The underlying goal is to compare the actual power decay with the
% inverse square law associated with a spherical wave in free field.
% 
% We do this comparison on active parts of the spectrum only, to
% eliminate the possible bias that would be introduced by bg noise.
% 
% This analysis requires a lapel WAV file, one or many distant
% microphone WAV files, and a ground-truth file, as for exemple in the
% recordings "seq01-1p-0000" and "seq03-1p-0000" of the AV16.3 Corpus
% ( http://mmm.idiap.ch/av16.3_v6 ).
% 
% Both input "p" and output "out" are structures.
% 
% The fields of "p" are the parameters of this function.
% Mandatory parameters:
% 
%   p.static_gt_filename        :   string
% 
%   p.lapelmic_filename         :   string
%
%   p.distantmic_filename_list  :   cell array of N_dist_mic strings
%
%   p.mic_list                  : integer array, 
%   the corresponding mic. indices in the "static_gt" structure.
%
%   p.framelength_sec           :   duration of frame, in seconds
%
% Optional parameters:
%   -> check the beginning of this code ( dbtype powerdecay_analysis ).
%
% The output structure "out" includes:
%
%   out.p : the input parameters, including the WAV filenames
%   specified by the user, for the lapel and the distant microphones
%   (see above). ".p" also includes optional parameters filled with
%   default values.
%
% and a structure array "out.loc", with as many elements as speaker
% locations. Each element has the following fields:
%
%   out.loc( i_loc ).lapelmic : a structure giving information about
%   the part of the lapel signal that we used, for this speaker
%   location.
%
%   out.loc( i_loc ).lapelmic.seg_sample_mat : a 2 x N_seg matrix of
%   integers (sample indices in the WAV file). Each column gives the
%   beginning and the end of each of the N_seg speech segments, for
%   this speaker location.
%
%   out.loc( i_loc ).lapelmic.w : a framelength x N_frames matrix of
%   weights, used to compute the narrowband statistics.
%
%   out.loc( i_loc ).lapelmic.wtf : a framelength x N_frames matrix of
%   weights, used to compute the fullband statistics.
%
%   out.loc( i_loc ).distantmic : a structure array of N_dist_mic
%   elements.
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).i_mic : integer, index
%   of the microphone in the microphone array.
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).wavfilename : string, 
%   name of the WAV file corresponding to "i_mic".
%
% <<< statistics computed on the whole spectrum (speech AND silence) >>>
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).mean_power_ratio_dB:
%   framelength x 1 vector, containing the narrowband average IN DB
%   DOMAIN of the power ratio ( P_distant / P_lapel ). This is equivalent
%   to a geometric mean in linear domain. NOTE THAT THE AVERAGE IS
%   COMPUTED ON ALL FRAMES (SPEECH AND SILENCE).
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).std_power_ratio_dB:
%   framelength x 1 vector of standard deviation values in dB domain,
%   similarly to ".mean_power_ratio_dB".
%
% <<< for each frame, the fullband power ratio >>
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).fullband_power_ratio_dB:
%   1 x N_frames vector of dB values: 10 * log10( P_distant / P_lapel ).
%
% <<< statistics computed ON SPEECH ONLY >>>
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).fullband_mean_power_ratio_dB:
%   scalar, average of ".fullband_power_ratio_dB", ON SPEECH FRAMES ONLY.
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).fullband_std_power_ratio_dB: 
%   scalar, standard deviation of ".fullband_power_ratio_dB", ON
%   SPEECH FRAMES ONLY.
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).w_mean_power_ratio_dB:
%   framelength x 1 vector, narrowband WEIGHTED average IN DB DOMAIN,
%   of the narrowband power ratio, where the weight is proportional to
%   1 on SPEECH (frame, frequency) and 0 on silent (frame, frequency).
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).w_std_power_ratio_dB:
%   framelength x 1 vector, narrowband WEIGHTED standard deviation, IN
%   DB DOMAIN, of the narrowband power ratio, where the weight is
%   proportional to 1 on SPEECH (frame, frequency) and 0 on silent
%   (frame, frequency).
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).wtf_mean_power_ratio_dB:
%   scalar, global weighted average, IN DB DOMAIN, of the narrowband
%   power ratio. The weight is proportional to 1 on SPEECH (frame,
%   frequency) and 0 on silent (frame, frequency).
%
%   out_loc( i_loc ).distantmic( i_dist_mic ).wtf_std_power_ratio_dB:
%   scalar, global weighted standard deviation, IN DB DOMAIN, of the
%   narrowband power ratio. The weight is proportional to 1 on SPEECH
%   (frame, frequency) and 0 on silent (frame, frequency).
%
% By Guillaume Lathoud, 2006  (email address in "my_email_address.jpg")
  
  if nargin < 1
    error( [ mfilename ' needs at least one input argument.' ] );
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 1 ) Get parameters from the user and the files
  
  % Check for mandatory input parameters

  required_param = { 'static_gt_filename', 'lapelmic_filename', 'distantmic_filename_list', 'mic_list', 'framelength_sec' };
  check_param( required_param, fieldnames( p ) );
  
  % Fill in default values for optional parameters

  p_default = [];

  p_default.frameshift_sec  = p.framelength_sec / 2;
  p_default.preemp          = 0.97;
  p_default.do_zero_mean    = 1;
  
  p_default.f_low  = 50;  % Value in Hz
  p_default.f_high = +Inf; % Value in Hz, or +Inf for Nyquist frequency
  p_default.n_narrowband = 80; % Number of narrowbands between "f_low" and "f_high"

  % Upsampling factor for the cross-correlation
  p_default.upfactor = 20;
  p_default.filterorder = [];

  % Tell me whether I should save the cross-correlograms
  % (this takes a lot of space)
  p_default.save_gcc = 1;
  
  % To get a reasonable starting point for delay search,
  % we assume a rough speed of sound in the air
  p_default.c = 342;  % in m/s

  % Assume that the lapel is 15 cm away from the mouth
  p_default.lapelmic_d = 0.15;
  
  % Proportion of the signals that we'll use.
  % We drop the queue of each utterance,
  % since it is usually weaker than the rest, 
  % (at least in average),
  % and therefore more sensitive to reverberations.
  p_default.utterance_prop = 0.75;

  % To limit the memory usage,
  % we limit the range of delays
  % around the theoretical free-field delay
  p_default.delay_min = -75;   % in sampling periods
  p_default.delay_max = +125;   % in sampling periods
  
  % Verbosity
  p_default.verbose = 1;

  % Additional parameters
  p_default.filterorder     = [];  % see below
  p_default.filterhalforder = [];  % see below
  p_default.B               = [];  % see below

  p = fill_default( p, p_default );
  
  %%%

  if isempty( p.filterorder )
    p.filterorder = 6 * p.upfactor;
  end

  if mod( p.filterorder, 2 )
    error( [ mfilename ': p.filterorder must be an even number.' ] );
  end

  if isempty( p.filterhalforder )
    p.filterhalforder = p.filterorder / 2;
  end

  if isempty( p.B )
    p.B = p.upfactor * fir1( p.filterorder, 1 / p.upfactor );
  end
  
  % Get the sampling frequency
  
  [ tmp, p.fs ] = wavread( p.lapelmic_filename, 'size' );
  
  p.f_high = min( p.fs / 2, p.f_high );
  
  % Parameters for the comparison
  
  p.framelength = round( p.framelength_sec * p.fs );
  p.frameshift  = round( p.frameshift_sec * p.fs );
  
  %%%

  if p.verbose
    disp( ' ' );
    disp( [ mfilename ' parameters:' ] );
    disp( p );
  end

  %%%

  % Bands in which we'll split the signal
  
  p.narrowbandwidth  = ( p.f_high - p.f_low ) / p.n_narrowband;
  p.narrowbandmat    = p.f_low + [ 0:p.n_narrowband-1; 1:p.n_narrowband ] * p.narrowbandwidth;
  p.narrowbandcenter = mean( p.narrowbandmat, 1 );
  
  % Create the corresponding filters
  
  p.narrowbandfilter = repmat( NaN, p.filterorder + 1, p.n_narrowband );

  for a = 1:p.n_narrowband
    
    if abs( p.narrowbandmat( 1, a ) ) < eps
      p.narrowbandfilter( :, a ) = fir1( p.filterorder, p.narrowbandmat( 2, a ).' / p.fs * 2, 'DC-1' ).';
    elseif abs( p.narrowbandmat( 2, a ) - p.fs / 2 ) < eps
      p.narrowbandfilter( :, a ) = fir1( p.filterorder, p.narrowbandmat( 1, a ).' / p.fs * 2, 'DC-0' ).';
    else
      p.narrowbandfilter( :, a ) = fir1( p.filterorder, p.narrowbandmat( :, a ).' / p.fs * 2, 'DC-0' ).';
    end    

  end

  % Move to the frequency domain
  
  p.narrowbandH  = fft( p.narrowbandfilter, 2 * p.framelength, 1 );
  p.narrowbandH2 = real( p.narrowbandH .* conj( p.narrowbandH ) );  % Real to remove possible imaginary epsilons
  
  p.narrowbandH2_GCC = p.narrowbandH2 > 1e-3;  % Necessary to filter GCC-PHAT
  
  % Based on the 3D ground-truth, estimate the relative distance between
  % lapel locations and microphones in the array.
  
  [ p.d, p.static_gt, p.p3d, p.g ] = estimate_dist( p.static_gt_filename );

  p.n_loc = size( p.d, 2 );

  p.d_diff = p.d - p.lapelmic_d;
  
  p.s_diff = p.d_diff / p.c * p.fs;  % time delay between lapel location and array microphones, in samples

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % ( 2 ) Start looping through speaker locations

  n_dist_mic = numel( p.mic_list );

  if p.verbose
    disp( ' ' );
    disp( sprintf( [ mfilename ': iterating through %d speaker locations.' ], p.n_loc ) );
    chrono = chrono_start( p.n_loc );
  end

  % Structure array where we'll store the output
  
  loc = [];
  
  for i_loc = 1:p.n_loc
    
    %%%%%%%%%%
    % ( 2.1 ) Define speech segments and get the lapel signal for this speaker location
    
    % Which speech segments
    
    i_seg = find( p.static_gt.pos_ind == i_loc );
    
    seg_sample_mat = zeros( 2, 0 );
    
    for qq = i_seg(:).'
      
      seg = p.static_gt.sp_seg( :, qq );
      seg = [ seg(1) ;  seg( 1 ) + p.utterance_prop * diff( seg ) ];
      seg_sample = round( seg(:) * p.fs );
      
      seg_sample_mat = [ seg_sample_mat  seg_sample ];
      
    end
    
    % Get the lapel signal
    
    lapelmic.seg_sample_mat = seg_sample_mat;
    lapelmic.S              = wavreadconcatspectrum( p.lapelmic_filename, seg_sample_mat, p );
    
    % Prepare some weights for computations below
    % ( it is meaningless to relate lapel and mic array in parts of the T-F plane where the speaker is silent )
    %
    % -> on speech only
    
    tmp      = abs( lapelmic.S );
    a_min    = min( tmp( find( tmp > 0 ) ) );
    if isempty( a_min )
      a_min = eps;
    end
    a_thr    = exp( mean( log( max( a_min, tmp ) ), 2 ) );
    tmp      = ( tmp > repmat( a_thr, 1, size( tmp, 2 ) ) );

    % To compute statistics for each band
    
    lapelmic.w = tmp ./ ( eps + repmat( sum( tmp, 2 ), 1, size( tmp, 2 ) ) );
    
    % To compute statistics for the whole T-F plane
    
    lapelmic.wtf = tmp / ( eps + sum( tmp(:) ) );
    
    %%%%%%%%%%
    % ( 2.2 ) Fullband cross-correlation analysis between the lapel and each microphone.
    %         Narrowband cross-correlation analysis between the lapel and each microphone.
    
    distantmic = [];
    
    for zzz = 1:n_dist_mic
      
      i_mic       = p.mic_list( zzz );
      wavfilename = p.distantmic_filename_list{ zzz };

      % Store some info

      newmic.i_mic       = i_mic;
      newmic.wavfilename = wavfilename;
      
      % Get the microphone array signal
      
      newmic.seg_sample_mat = lapelmic.seg_sample_mat + round( p.s_diff( i_mic, i_loc ) );
      newmic.S              = wavreadconcatspectrum( wavfilename, newmic.seg_sample_mat, p );
      
      %%% First, raw ratios for each frame and each (frame, narrowband)

      % Compare the spectra in a bandwise manner
      
      newmic.mag_ratio      = abs( newmic.S ) ./ ( eps + abs( lapelmic.S ) );
      newmic.power_ratio_dB = 20 * log10( newmic.mag_ratio );
      
      newmic.mean_power_ratio_dB = mean( newmic.power_ratio_dB, 2 );
      newmic.std_power_ratio_dB  = std( newmic.power_ratio_dB, [], 2 );
      
      % Also compare the fullband powers
      
      tmp_1 = sum( abs( newmic.S ) .^ 2, 1 );
      tmp_2 = sum( abs( lapelmic.S ) .^ 2, 1 );
      
      newmic.fullband_power_ratio_dB = 10 * log10( tmp_1 ./ ( eps + tmp_2 ) );
      
      %%% Second, average statistics using weighting
      %%% --> we compare powers ONLY on parts of the T-F plane where the speaker is active
      
      % Fullband ratio: only consider speech frames
      
      a_thr = exp( mean( log( tmp_1( find( tmp_1 > 0 ) ) ) ) );
      ind   = find( tmp_1 > a_thr );
      
      newmic.fullband_mean_power_ratio_dB = mean( newmic.fullband_power_ratio_dB( ind ), 2 );
      newmic.fullband_std_power_ratio_dB  = std( newmic.fullband_power_ratio_dB( ind ), [], 2 );
      
      % Narrowband ratio: give more weight to (band,time) with large magnitude
      % ( it is meaningless to relate lapel and mic array in parts of the T-F plane where the speaker is silent )
      % ( Above, we have already made sure that weights sum to one ).
      
      newmic.w_mean_power_ratio_dB = sum( lapelmic.w .* newmic.power_ratio_dB, 2 );
      
      tmp = repmat( newmic.w_mean_power_ratio_dB, 1, size( newmic.power_ratio_dB, 2 ) );
      newmic.w_std_power_ratio_dB  = sqrt( sum( lapelmic.w .* ( ( newmic.power_ratio_dB - tmp ) .^ 2 ), 2 ) );
      
      % And finally a global SNR measure across the whole T-F plane
      
      newmic.wtf_mean_power_ratio_dB = sum( lapelmic.wtf(:) .* newmic.power_ratio_dB(:) );
      newmic.wtf_std_power_ratio_dB  = sqrt( sum( lapelmic.wtf(:) .* ( ( newmic.power_ratio_dB(:)- newmic.wtf_mean_power_ratio_dB ).^2 ) ) );
      
      % Remove the mic spectrum, it adds up to too much space overall
      
      newmic.S = [];
      newmic.mag_ratio = [];
      newmic.power_ratio_dB = [];

      %%% Store the results
            
      distantmic = cat( 1, distantmic, newmic );
      
    end  % loop through microphones

    % Remove spectra

    distantmic = rmfield( distantmic, 'S' );
    distantmic = rmfield( distantmic, 'mag_ratio' );
    distantmic = rmfield( distantmic, 'power_ratio_dB' );
    lapelmic   = rmfield( lapelmic, 'S' );
    
    % Save the result
    
    loc( i_loc ).lapelmic   = lapelmic;
    loc( i_loc ).distantmic = distantmic;

    % Verbosity
    
    if p.verbose
      chrono = chrono_check( chrono, i_loc );
    end
    
  end   % Of loop through speaker locations

  if p.verbose
    chrono_stop( chrono );
  end

  out = [];
  out.p = p;
  out.loc = loc;
  
  return;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ d, static_gt, p3d, g ] = estimate_dist( static_gt_filename )

  % Load spatial and temporal ground-truth
  % ( one speaker at several locations )
  
  load( static_gt_filename, 'static_gt' );
  
  % Position of the lapel w.r.t. middle of the two microphone arrays
  
  p3d = static_gt.p3d;
  
  % Geometry of each array (copied from "seq_av163_modified.m" in the AV16.3 directory)

  a = (4:11)*pi/4;
  g = 0.1* [cos(a);sin(a);0*a];

  g1 = g; g1( 2,: ) = g1( 2,: ) + 0.4;
  g2 = g; g2( 2,: ) = g2( 2,: ) - 0.4;

  g = [ g1  g2 ];
  
  % Estimate the distance between each microphone in the array and each lapel location

  n_loc = size( p3d, 2 );
  n_mic = size( g, 2 );

  d = repmat( NaN, n_mic, n_loc );

  for a = 1:n_mic
    
    delta = p3d( 1:3, : ) - repmat( g( :, a ), 1, n_loc );

    d( a, : ) = sqrt( sum( delta .^ 2, 1 ) );

  end
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ S, fs, nbits ] = wavreadconcatspectrum( lapelmic_filename, seg_sample_mat, p )

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

  S = zeros( 2 * p.framelength, 0 );
  
  for a = 1:size( seg_sample_mat, 2 )

    [ x, fs, nbits ] = wavread( lapelmic_filename, seg_sample_mat( :, a ).' );

    nx = numel( x );
    
    % Full frames only
    new_nx = p.framelength + floor( ( nx - p.framelength ) / p.frameshift ) * p.frameshift;

    x = x( 1:new_nx );

    S = [ S  my_spectrum( x, p.framelength, p.frameshift, p.preemp, p.do_zero_mean ) ];
  
  end
  