function X = wavelet(x,T,a,wwin)
%
% WAVELET - performs a wavelet transform
%
%    The command X = wavelet(x,T,a) will return the wavelet transform
%  of the signal contained in the single vector x.  The returned
%  argument X is a two dimmensional function X(T,a) where T is time
%  delay and a is time scale which can be roughly thought of as 
%  frequency.  X(T,a) is best visualized with the MESH command.  The
%  wavelet function used is a common wavelet function defined in
%  the file WAVEWIN.M.  The variable a is a vector containing
%  the desired range of time scale.  The variable T is a vector of
%  INTEGERS spanning the desired time delay range.  Since the signal
%  in x is discrete, no more time resolution can be obtained.  Scale
%  a, however may be of any range.  Try using logslpace to create
%  useful a vectors.   
%
%    If the optional string argument 'window' is specified as in
%  y=wavelet(x,T,a,'window') then the wavelet function contained in
%  the .m-file 'window' is used instead of WAVEWIN. 
%
%  NOTE: When writing custom wavelet windowing functions, it is
%        important that they conform to certain mathematical 
%        criteria.  You should be familiar with these.  You may wish
%        to use WAVEWIN.M as a template for your wavelet.
%
%  The wavelet transform is defined as follows:
%
%                              +inf
%                       1     |^   
%             X(T,a) = ----  |  x(t)ˇgamma[(t-T)/a]ˇdt
%                      a^˝  V
%                            -inf
%
%               where gamma(t) is a suitable wavelet function.
%
x=x(:)';
T=T(:)';
a=a(:)';
n=length(x);

if nargin == 3,
  wwin='wavewin';
end

%                Pad appropriately with zeros
%                --- ------------- ---- -----
if T(1) < 0
  pad1 = zeros(1,-T(1));
end
if T(length(T)) > (n-1)
  pad2 = zeros(1,length(T)-n+1);
end

x = [pad1 x pad2];
 

%              Perform the Wavelet Transform
%              ------- --- ------- ---------
X=zeros(1,length(x));
Xtemp=[];
for i = 1:length(a),
 for j = 1:length(T),
   y=feval(wwin,x,T(j),a(i));
   Xtemp=cumsum(y)/sqrt(a(i));
 end
 X = [X;Xtemp];
 Xtemp=[];
end

%               Now Plot It
%               --- ---- --
surf(X),shading('interp'),view([0 90])
title('Wavelet Transform Scalogram'),xlabel('Time'),ylabel('Scale')