awt.sci

小波分解源代码

function awt = AWT(x,nvoice,Name,sigma2,par2)
// AWT -- Analytical Wavelet Transform
//  Usage
//    awt = AWT(x,nvoice,Name,sigma2,par2);
//  Inputs
//    x	     signal, could be complex
//    nvoice  number of voices
//            default = 12
//    Name    Type of Window function, so far only "Gaussian" window available
//    sigma2    first parameter, for Gaussian window, it could be variance,
// 	     default = 1
//   par2    another parameter, for Gaussian, it's xi(shifting frequency)
//            default = 5
//  Outputs
//    awt     Analytical wavelet transform of x
//  Side Effects
//  Description
//    Coarsest scale is '2', finest scale is (log2(n)-5);
//    awt is a n by nscale matrix;
//  Algorithm
//     
//  Examples
//    n = 1024;
//    t = ((1:n)./n - .5).*%pi;
//    f = exp(%i*t.^2);//x=2*sin(2*t)+3*cos(5*t).*sin(2*t);
//    
//    awt = AWT(f);
//   
//    ix = [((n/2+1):n),(1:n/2)];
//   plot(abs(awt(ix,1)))	// should looks Gaussian
//    mtlb_hold on
//    plot(abs(awt(ix,13)))	// should looks Gaussian
//    plot(abs(awt(ix,25)))	// should looks Gaussian
//  See Also
//    ImageCWT
//  References
//    Mallat, "A Wavelet Tour of Signal Processing";
//             4.4.3 Analytical Wavelet Transform
//
//  Copyright Aldo I Maalouf

  n = length2(x);
[lhs,rhs]=argn();
// set default values
  if  rhs < 4,
	sigma2 = 1; xi = 5;
  end;
 if  rhs < 3,
	Name = 'Gabor';
  end;
 if  rhs < 2,
 	nvoice = 12;
  end;

// Fourier Transform of signal
  dftx = mtlb_fft(x(:));
  
  omega   = [ (0: (n/2)) (((-n/2)+1):-1) ].* (2*%pi/n);

  noctave = floor(log2(n))-4;
  nscale  = nvoice .* noctave;

  awt = zeros(n,nscale);

  kscale  = 1;
  scale   = 2;

  for jo = 1:noctave,
     for jv = 1:nvoice,
        s =  scale .* (2^(jv/nvoice));
        if string(Name)=='Gabor';
        
      	      fre =  (s .* omega  - xi);
              fre = fre(:);
	      omega = omega(:);
	      Psi = realpow(4.*%pi.*sigma2,1/4)*exp(-sigma2/2*fre.*fre);
              Psi = Psi .* (omega>0);
        end; 
	awt(1:n,kscale) = mtlb_ifft(dftx.*Psi);   
	
	kscale = kscale + 1;
     end;
     scale  = scale .*2;
  end;

endfunction