function new_x = overlap_add2( X, filter_mat, framelength, frameshift, preemp )

% function new_x = overlap_add( X, filter_mat, framelength, frameshift, preemp )
%
% overlap-add reconstruction of a signal, from its COMPLEX-domain spectrogram "X".
%
% Inputs:
%
%   X (mandatory): 2*framelength by nframes matrix of complex values,
%     symmetric spectrogram (as given e.g. by a FFT).
%
%   filter_mat (optional): 2*framelength by nframes matrix of real or
%     complex values.  It defines a filter that you want to apply on
%     X.  For example you may give a matrix of ones and zeroes, if you
%     only want to keep certain (time, frequency) points of your
%     spectrogram.
%
%    framelength (mandatory): integer, length of your time frame, in samples.
%
%    frameshift (mandatory): integer, shift between two time frames, in samples.
%
%    preemp (optional, default 0.97): pre-emphasis coefficient.
%
%
% Output:
%
%    new_x: vector of real values, time-domain waveform reconstructed
%      from the (optionally filtered) complex spectrogram X.
%
% By Guillaume Lathoud 2005 - lathoud@idiap.ch
% 
% Thanks to Oliver Masson for all the explanations.


if nargin < 4
   error( [ mfilename ' needs at least 4 input arguments.' ] ); 
end

if ~exist( 'preemp', 'var' )
    preemp = 0.97;
end

% XXX slightly better this way
% (on tests we found it slightly closer to the original signal)
if isempty( filter_mat )
  filter_mat = ones( size( X ) );
end

if ~isempty( filter_mat )
  if numel( unique( [ size( X, 1 ) size( filter_mat, 1 ) 2*framelength ] ) ) > 1
    error( [ mfilename ': inconsistent frame length.' ] ); 
  end
  
  nframes = unique( [ size( X, 2 ) size( filter_mat, 2 ) ] );
  if numel( nframes ) > 1
    error( [ mfilename ': inconsistent number of frames.' ] );     
  end 
else
  if numel( unique( [ size( X, 1 ) 2*framelength ] ) ) > 1
    error( [ mfilename ': inconsistent frame length.' ] ); 
  end  
  nframes = size( X,2 );
end

% overlap-add construction of a signal

% scaling + masking

if ~isempty( filter_mat )

  new_filter = ifft( filter_mat, 2*framelength, 1 );
  new_filter = new_filter( [ end-framelength/2+1:end 1:framelength/2 ], : );
  new_filter = diag( sparse( blackman( framelength, 'periodic' ) ) ) * new_filter;
  new_filter( framelength+1:2*framelength, : ) = 0;
  
  % Apply the filter in a correct manner, in the frequency domain
  new_X = X .* fft( new_filter, 2*framelength, 1 );
  
else
  
  new_X = X;

end

% overlap-add
nsamples = framelength + ( nframes-1 ) * frameshift;
new_x = zeros( nsamples+4*framelength, 1 );
the_long_win = hamming( 2*framelength, 'periodic' );
for a = 1:nframes
      
    % The "real" is for very very small numerical imprecision
    % in order to remove some 1e-10-like imaginary part
    a_frame = real( ifft( new_X( :,a ), 2 * framelength )) ./ the_long_win;
    start_sample = 1 + (a-1) * frameshift;
    end_sample   = start_sample + 2 * framelength - 1;
    ind = start_sample:end_sample;
    new_x( ind ) = new_x( ind ) + a_frame;
end

% Finally remove pre-emphasis
new_x = filter( 1, [ 1 -preemp], new_x );
  
if ~isempty( filter_mat )
  % Compensate the filter's delay: cut the crap at the beginning
  new_x = new_x( framelength/2+1:end );
end

% Cut off the crap at the end
new_x = new_x( 1:nsamples );

% Return a column vector
new_x = new_x(:);