function [x,t] = tfristft(tfr,t,h,trace);
%TFRISTFT Inverse Short time Fourier transform.
%	[X,T]=TFRSTFT(tfr,T,H,TRACE) computes the inverse short-time 
%	Fourier transform of a discrete-time signal X. This function
%	may be used for time-frequency synthesis of signals.
% 
%	X     : signal.
%	T     : time instant(s)          (default : 1:length(X)).
%	H     : frequency smoothing window, H being normalized so as to
%	        be  of unit energy.      (default : Hamming(N/4)). 
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR   : time-frequency decomposition (complex values). The
%	        frequency axis is graduated from -0.5 to 0.5.
%
%	Example :
%	 t=200+(-128:127); sig=[fmconst(200,0.2);fmconst(200,0.4)]; 
%        h=hamming(57); tfr=tfrstft(sig,t,256,h,1);
%        sigsyn=tfristft(tfr,t,h,1);
%        plot(t,abs(sigsyn-sig(t)))
% 
%	See also all the time-frequency representations listed in
%	the file CONTENTS (TFR*)

%	F. Auger, November 1996.
%
%	------------------- CONFIDENTIAL PROGRAM -------------------- 
%	This program can not be used without the authorization of its
%	author(s). For any comment or bug report, please send e-mail to 
% 	f.auger@ieee.org

if (nargin<3),
 error('At least 3 parameters required');
elseif (nargin==3),
 trace=0;
end;

[N,NbPoints]=size(tfr);
[trow,tcol] =size(t);
[hrow,hcol] =size(h); Lh=(hrow-1)/2; 

if (trow~=1),
 error('T must only have one row'); 
elseif (hcol~=1)|(rem(hrow,2)==0),
 error('H must be a smoothing window with odd length');
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
elseif (tcol~=NbPoints)
 error('tfr should have as many columns as t has rows.');
end; 

Deltat=t(2:tcol)-t(1:tcol-1); 
Mini=min(Deltat); Maxi=max(Deltat);
if (Mini~=1) & (Maxi~=1),
 error('The tfr must be computed at each time sample.');
end;

h=h/norm(h);

if trace, disp('Inverse Short-time Fourier transform'); end;
tfr=ifft(tfr);

x=zeros(tcol,1);

for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 valuestj=max([1,icol-N/2,icol-Lh]):min([tcol,icol+N/2,icol+Lh]);
 for tj=valuestj,
  tau=icol-tj; indices= rem(N+tau,N)+1; 
  % fprintf('%g %g %g\n',tj,tau,indices);
  x(icol,1)=x(icol,1)+tfr(indices,tj)*h(Lh+1+tau);
 end;
 x(icol,1)=x(icol,1)/sum(abs(h(Lh+1+icol-valuestj)).^2);
end;

if trace, fprintf('\n'); end;