ctwt.m

这是伯克里wavelet transforms一书中的例子的代码

function c = ctwt(sig,ts,ain,wopt)
%Routine to compute the scalogram of the given sequence.
%
%c = ctwt(seq,ts,a,opt) computes the scalogram of the input sequence,seq
%that has been sampled at the rate, ts. The scalogram is computed for the
%range, -a to a with a step size of 0.1. The variable, opt chooses the  
%wavelet to be used. The contour plot of the scalogram is displayed and
%the array containing the scalogram is returned. 
%
%The variables are:
%seq      : The input sequence.
%ts       : Sampling period of the input sequence.
%a        : The largest scale. Scalogram is computed from -a to a in
%           steps of 0.1
%opt      : This chooses the wavelet. The options are as follows:            
%           1 - Haar wavelet
%           2 - Mexican Hat wavelet
%           3 (or >3) - Morlet wavelet.
%
%Example:
%
%c = ctwt(data,0.1,2,3)
%
%Will compute the scalogram of the given one-dimensional signal, data, 
%assuming that it has been sampled every 0.1 units. The scalogram will be 
%computed for the range of scales -2 to 2 in steps of 0.1 using the Morlet
%wavelet (wavelet option 3). The resulting scalogram will be returned in the 
%array, c.
%
%Chapter 1 gives more information on the continuous wavelet transform and
%the scalogram.
%
%Author: Ajit S. Bopardikar
%Copyright (c) 1998 by Addison Wesley Longman, Inc.
%

[r,c] = size(sig);
if (r ~=1)
   sig = sig';
end %endif
 
l =length(sig);
ain = round(ain*10)/10;

if (wopt == 1) %Haar wavelet

  c = []; %Initialize

  for(a=abs(ain):-0.1:-abs(ain))

     if(a > 0)
       ha = haar(a,ts);
       lh = length(ha);
       ct = -conv(sig,ha);
        c = [c; ct(lh:lh+l-1)];
     elseif (a < 0)
       ha = haar(a,ts);
       lh = length(ha);
       ha = [zeros([1 lh]) ha]; 
       ct = -conv(sig,ha);
        c = [c; ct(lh:lh+l-1)];
     elseif (a == 0)
       [rows,cols]=size(c);
       c = [zeros([1 cols]); c];
     end %endif 
  end  %endfor
     
else  %Mexican hat wavelet ot Morlet wavelet
  cker = []; %initialize the kernel... 

  for a=(abs(ain):-0.1:-abs(ain))  %This procedure is for the options of Morlet or 
                       % Mexican hat wavelet.
     if (a~=0)
       if(wopt == 2)            %Mexican hat wavelet
         f = mexhat(a,ts,5*ain);
       elseif (wopt >= 3)       %Morlet wavelet
         f = morlet(a,ts,5*ain);
       end %endif 
       cker = [f;cker];
    else  %a=0
        [rows,cols] = size(cker);
        cker = [zeros([1 cols]);cker];
    end %endif       
  end %endfor

  c = conv2(sig,cker);  %Compute the CTWT...
  c = c(:,round(cols/2):round(cols/2)+l);
end %endif

[r,co] = size(c);
y =abs(ain):-0.1:-abs(ain);
x = (0:co-1)*ts;
c=abs(c).^2;
figure;contour(x,y,c,20);
xlabel('b');ylabel('a');