function [tfr,t,f]=tfrsp(x,t,N,h,tr)
tfr=[];f=[];
[nargout,nargin] = argn(0)
//TFRSP Spectrogram time-frequency distribution.
// [TFR,T,F]=TFRSP(X,T,N,H,TRACE) computes the Spectrogram 
// distribution of a discrete-time signal X. 
// 
// X     : signal.
// T     : time instant(s)          (default : 1:length(X)).
// N     : number of frequency bins (default : length(X)).
// H     : analysis 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 representation. When called without 
//         output arguments, TFRSP runs TFRQVIEW.
// F     : vector of normalized frequencies.
// 
// Example :
//  sig=fmlin(128,0.1,0.4);
//  h=tftb_window(17,'Kaiser'); tfrsp(sig,1:128,64,h,1);
// 
//        [tfr,t,freq]=tfrsp(sig,1:128,64,h,1); plot(fftshift(freq),fftshift(tfr(:,100)))
// 
// See also all the time-frequency representations listed in
//  the file CONTENTS (TFR*)
 
// F. Auger, May-August 1994, July 1995.
// Copyright (c) 1996 by CNRS (France).
// 
// ------------------- 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.rog 

if nargin==0 then
  error('At least 1 parameter required');
end
[xrow,xcol] = size(x);

if xcol>1 then
  error('Input signal X must be a column vector');
end
 
if nargin<1 then
  error('At least 1 parameter is required');
elseif nargin<=2 then
  N = xrow;
end
 
hlength = floor(N/4);
hlength = hlength+1-(hlength-fix(hlength/2)*2);
 
if nargin==1 then
  t = 1:xrow;
  h = window('hn',hlength)';
  tr = 0;
elseif nargin==2|nargin==3 then
  h = window('hn',hlength)';
  tr = 0;
elseif nargin==4 then
  tr = 0;
end
 
if N<0 then
  error('N must be greater than zero');
end
[trow,tcol] = size(t);
if xcol==0|xcol>2 then
  error('X must have one or two columns');
elseif trow~=1 then
  error('T must only have one row');
elseif (2^ceil(log(N)/log(2)))~=N then
  printf('For a faster computation, N should be a power of two\n');
end
 
[hrow,hcol] = size(h);
Lh = (hrow-1)/2;

if hcol~=1|(hrow-fix(hrow./2).*2)==0 then
  error('H must be a smoothing window with odd length');
end
 
tfr = zeros(N,tcol);
if tr, printf('Spectrogram'); end;
for icol = 1:tcol
  ti = t(icol);
  tau = -min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
  dummy = N+tau
  indices = dummy-fix(dummy./N).*N+1;
  if tr, disprog(icol,tcol,10); end;
  tfr(indices,icol) = x(ti+tau) .* conj(h(Lh+1+tau))/norm(h(Lh+1+tau));
  tfr(:,icol)=fft(tfr(:,icol),-1);
end
 
tfr = abs(tfr).^2;
if tr, printf('\n'); end;
 
if nargout==3 then
  if (N-fix(N./2).*2)==0 then
    f = [0:N/2-1,-N/2:-1]'/N;
  else
    f = [0:(N-1)/2,-(N-1)/2:-1]'/N;
  end
end
 
