function show_powerdecay_summary( powerdecay_result_filename, loc_ind, mic_ind )

% Show the result of "powerdecay_analysis_wrapper.m"
% in a graphical way:
% 
%   A summary picture with the mean and standard deviation
%   of the power ratio, as a function of the speaker distance.
%   
%   y = power ratio between distant mic. and lapel mic., in dB
%   x = difference of distance between lapel and distant microphone
%
% Refer to "powerdecay_analysis_wrapper.m" and "powerdecay_analysis.m"
% for more details.
%
% By Guillaume Lathoud, 2006  (email address in "my_email_address.jpg")

  if nargin < 1
    error( [ mfilename ' needs 1 input argument.' ] );
  end
  
  if ~exist( 'loc_ind', 'var' )
    loc_ind = [];
  end

  if ~exist( 'mic_ind', 'var' )
    mic_ind = [];
  end

  % load the result
  disp( [ mfilename ' is loading "' powerdecay_result_filename '"...' ] );
  load( powerdecay_result_filename, 'out' );
  
  % Accumulate "x", "y1" and "y2" data, as defined above
  x     = [];
  y1     = [];
  y1_bar = [];
  y2     = [];
  y2_bar = [];

  n_loc = numel( out.loc );
  
  disp( sprintf( [ mfilename ': iterating through the %d speaker locations...' ], n_loc ) );
  
  if isempty( loc_ind )
    i_loc_list = 1:n_loc;
  else
    i_loc_list = intersect( 1:n_loc, loc_ind );
  end
  
  for i_loc = i_loc_list
    
    n_mic = numel( out.loc( i_loc ).distantmic );
    
    if isempty( mic_ind )
      i_mic_list = 1:n_mic;
    else
      i_mic_list = intersect( 1:n_mic, mic_ind );
    end
    
    for i_mic = i_mic_list
      
      % x: difference of speaker distance between lapel and distant microphone, in meters
      
      x     = [ x     out.p.d_diff( out.loc( i_loc ).distantmic( i_mic ).i_mic, i_loc ) ];
      
      % y1: RMS of normalized narrowband power ratio (on active bands only)
      
      y1     = [ y1     out.loc( i_loc ).distantmic( i_mic ).wtf_mean_power_ratio_dB ]; 
      
      y1_bar = [ y1_bar out.loc( i_loc ).distantmic( i_mic ).wtf_std_power_ratio_dB ]; 
      
      % y2: RMS of normalized narrowband power ratio (on active bands only)
      
      y2     = [ y2     out.loc( i_loc ).distantmic( i_mic ).fullband_mean_power_ratio_dB ]; 
      
      y2_bar = [ y2_bar out.loc( i_loc ).distantmic( i_mic ).fullband_std_power_ratio_dB ]; 
      
    end

  end

  %%%

  % Sort the data by increasing speaker distance

  [ x, ind ] = sort( x );
  
  y1     = y1( ind );
  y1_bar = y1_bar( ind );

  y2     = y2( ind );
  y2_bar = y2_bar( ind );
  
  %%%

  [ hfig hplot ] = my_figure( x, y1, y1_bar, [ out.recording_name ' - narrowband' ] );

  % Graphical output to files
  [ p, n, e ] = fileparts( powerdecay_result_filename );
  
  % Dump eps to file (b/w)
  figure( hfig ); colormap( 0.1 + 0.9 * ( colormap( 'gray' ) .^ 3 ) );
  set( hplot, 'color', 'k' );
  out_filename = fullfile( p, [ n '-powerratio_narrowband-bw.eps' ] );
  print( hfig, '-depsc2', out_filename );
  disp( [ mfilename ': wrote out "' out_filename '".' ] );
  
  % Dump eps to file (color)
  figure( hfig ); colormap( jet );
  set( hplot, 'color', 'b' );
  out_filename = fullfile( p, [ n '-powerratio_narrowband-color.eps' ] );
  print( hfig, '-depsc2', out_filename );
  disp( [ mfilename ': wrote out "' out_filename '".' ] );

  %%%

  [ hfig hplot ] = my_figure( x, y2, y2_bar, [ out.recording_name ' - fullband' ] );

  % Graphical output to files
  [ p, n, e ] = fileparts( powerdecay_result_filename );
  
  % Dump eps to file (b/w)
  figure( hfig ); colormap( 0.1 + 0.9 * ( colormap( 'gray' ) .^ 3 ) );
  set( hplot, 'color', 'k' );
  out_filename = fullfile( p, [ n '-powerratio_fullband-bw.eps' ] );
  print( hfig, '-depsc2', out_filename );
  disp( [ mfilename ': wrote out "' out_filename '".' ] );
  
  % Dump eps to file (color)
  figure( hfig ); colormap( jet );
  set( hplot, 'color', 'b' );
  out_filename = fullfile( p, [ n '-powerratio_fullband-color.eps' ] );
  print( hfig, '-depsc2', out_filename );
  disp( [ mfilename ': wrote out "' out_filename '".' ] );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [ hfig, hplot ] = my_figure( x, y, y_bar, a_title )

  if nargin < 4
    error( [ mfilename ' needs 4 input arguments.' ] );
  end
  
  % Estimate what the inverse and inverse square laws would give
  
  % inverse law (no physical justification, just for comparison)
  
  d1_y_prop  = -10 * log10( x );   
  d1_y_const = mean( y - d1_y_prop );
  d1_y       = d1_y_prop + d1_y_const;
  
  % inverse square law (sperical wave propagating in a free field,
  % from a point source)
  
  d2_y_prop  = -20 * log10( x );   
  d2_y_const = mean( y - d2_y_prop );
  d2_y       = d2_y_prop + d2_y_const;
  
  % now try to "fit" the exponent of the law, in a MSE sense in dB domain.
  %
  % We start from the inverse square law and run Scaled Conjugate Gradient descent (Moller 93).
  
  global data
  data.x = x;
  data.y = y;
  
  ws = [];
  ws.w_1 = [ 2; d2_y_const ];
  ws.dE_E_fun = @my_cost_function;
  ws.cv_thr   = 1e-7;

  ws = run_scg_fast( ws );
  
  fit_alpha = ws.w_kp1( 1 );
  fit_const = ws.w_kp1( 2 );

  dfit_y  = -fit_alpha * 10 * log10( x ) + fit_const;
  
  % Now start plotting
  
  hfig = figure; hold on;
  
  % Second layer (on top): measurements
  
  hplot = errorbar( x, y, y_bar, 'k.' ); set( hplot, 'markersize', 20 ); 
  
  set( gca, 'FontSize', 20 );
  
  % First layer (below): d^{-1} and d^{-2} laws
  
  v = [ min( x )  max( x ) ];
  
  a_linewidth = 2;

%  hd1 = plot( x, d1_y, 'k-.', 'linewidth', a_linewidth );
  
  hd2 = plot( x, d2_y, 'm--', 'linewidth', a_linewidth );
  
  hdfit = plot( x, dfit_y, 'r-', 'linewidth', a_linewidth );

  % Annotate the graph
  
  title( strrep( a_title, '_', '\_' ) );
  
  xlabel( 'difference of distance (meters)' );
  ylabel( 'power ratio ( P_{distant} / P_{lapel} ) in dB' );
  
%  legend( [ hplot( [ 2 1 ] ).'  hd1  hd2 hdfit ], 'mean', 'std dev', '||\cdot||^{-1} law', '||\cdot||^{-2} law', [ '||\cdot||^{-' sprintf( '%f', fit_alpha ) '} law' ] );
  legend( [ hplot( [ 2 1 ] ).'   hd2 hdfit ], 'mean', 'std dev', '||\cdot||^{-2} law', [ '||\cdot||^{-' sprintf( '%f', fit_alpha ) '} law' ] );
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ dE, E ] = my_cost_function( w )
  
% function [ dE, E ] = my_cost_function( w )
%
% The data sits in "ws.data_x" and "ws.data_y".
% We are trying to find a function f( y ) = - alpha * 10 * log10( x ) + const
% that will best approximate the data in a mean square sense.
% 
% The present Matlab function evaluate the gradient "dE" of the cost "E",
% where E = mean( ( ws.data_y - f( ws.data_x ) .^ 2 )
  
  global data

  if nargin < 1
    error( [ mfilename '.my_cost_function() needs 1 input arguments.' ] );
  end
  
  % Extract the parameters

  if numel( w ) ~= 2
    error( [ mfilename '.my_cost_function(): input argument "w" must have 2 elements.' ] );
  end
  
  alpha = w( 1 );
  const = w( 2 );

  % Common term in both gradient and cost
  
  common_term = data.y + alpha * 10 * log10( data.x ) - const;

  % Compute the gradient
  
  dE = [ NaN; NaN ];
  
  dE( 1 ) = mean( 2 * 10 * log10( data.x ) .* common_term );  % with respect to "alpha"
  
  dE( 2 ) = mean( 2 *        ( -1 )        .* common_term );  % with respect to "const"
  
  if nargout > 1
    
    % Compute the cost
    
    E = mean( common_term .^ 2 );
    
  end
  
