function y = tffilter(tfr,x,t,trace);
%TFFILTER Time frequency filtering of a signal.
%       Y=TFFILTER(TFR,X,T,TRACE) filters the signal X
%	with a non stationary filter. 
%	
% 
%	X     : input signal (must be analytic).
%	T     : time instant(s)          (default : 1:length(X)).
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 (default : 0).
%	TFR   : Wigner-Ville distribution of the filter
%	        frequency axis is graduated from 0.0 to 0.5.
%
%	Example :
%	 Nt=128; t=1:Nt; sig=fmlin(Nt,0.05,0.3)+fmlin(Nt,0.2,0.45); 
%        sig(Nt/2)=sig(Nt/2)+8; figure(1);tfrwv(sig,t);
%        Nf=128;freqs=0.5*(0:Nf-1).'/Nf; 
%        for tloop=1:Nt, 
%         rate=0.2*(tloop-1)/(Nt-1);
%         H(:,tloop)=(0+rate<freqs).*(freqs<0.1+rate); 
%        end;
%        y=tffilter(H,sig,t,1);figure(2); tfrwv(y,t);
%       Example 2:
%	 Nt=128; t=1:Nt; sig=atoms(128,[64 0.25 sqrt(128) 1]);
%        figure(1);tfrwv(sig,t);
%        Nf=64;
%	 H=zeros(Nf,Nt);H(Nf/4+(-15:15),Nt/2+(-15:15))=ones(31);
%        y=tffilter(H,sig,t,1);figure(2); tfrwv(y,t);
%
% 
%	See also: tfristft

%	F. Auger, Jan 1997.
%
%	------------------- 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 
%       auger@ge44.univ-nantes.fr 

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

[N,NbPoints]=size(tfr);
[trow,tcol] =size(t);
[xrow,xcol] =size(x);

if (trow~=1),
 error('T must only have one row'); 
elseif (xcol~=1),
 error('X must only have one column'); 
elseif (2^nextpow2(N)~=N),
 fprintf('For a faster computation, N should be a power of two\n');
elseif (xrow~=NbPoints)
 error('tfr should have as many columns as X has rows.');
elseif (min(t)<1)|(max(t)>NbPoints),
 error('the values of T must be between 1 and xrow.');
end; 

if trace, disp('non stationary signal filtering'); end;

tfr=ifft(tfr);
y=zeros(tcol,1);

for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 ti=t(icol); 
 valuestj=max([1 2-ti ti-N/2+1]):min([NbPoints 2*NbPoints-ti ti+N/2]);
 %fprintf('%f %f %f\n',ti, min(valuestj),max(valuestj));

 for tj=valuestj,
  tmid = fix(0.5*(ti+tj)); 
  tdiff= fix(ti-tj); indice= rem(N+tdiff,N)+1;
  %fprintf('%g %g %g\n',tj,tmid,indice);
  y(icol,1)=y(icol,1)+tfr(indice,tmid)*x(tj);
 end;
end;

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