cwt.sci

小波分解源代码

function cwt = CWT(x,nvoice,wavelet,oct,scale)
// CWT -- Continuous Wavelet Transform
//  Usage
//    cwt = CWT(x,nvoice,wavelet,oct,scale)
//  Inputs
//    x        signal, dyadic length n=2^J, real-valued
//    nvoice   number of voices/octave
//    wavelet  string 'Gauss', 'DerGauss','Sombrero', 'Morlet'
//    octave   Default=2
//    scale    Default=4
//   Outputs
//   cwt      matrix n by nscale where
//             nscale = nvoice .* noctave
//  
//  Copyright Aldo I Maalouf
        [lhs,rhs]=argn();
	if rhs<4,
		oct = 2;
		scale = 4;
	end	
// preparation
	x = ShapeAsRow(x);
	n = length2(x);
	xhat = mtlb_fft(x);
	xi   = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*%pi/n);
	
// root
	omega0 = 5;
	

	noctave = floor(log2(n))-oct;
	nscale  = nvoice .* noctave;
	
	cwt = zeros(n,nscale);

	kscale  = 1;


	for jo = 1:noctave,
	    for jv = 1:nvoice,
		   qscale = scale .* (2^(jv/nvoice));
		   omega =  n .* xi ./ qscale ;
		   
		   if string(wavelet)=='Gauss',
				wind = exp(-omega.^2 ./2);
		
		elseif string(wavelet)=='DerGauss',
                                wind = %i.*omega.*exp(-omega.^2 ./2);
		   
		elseif string(wavelet)=='Sombrero',	
				wind = (omega.^2) .* exp(-omega.^2 ./2);
		   
		elseif string(wavelet)=='Morlet';	
				wind = exp(-(omega - omega0).^2 ./2) - exp(-(omega.^2 + omega0.^2)/2);
		   end
		   // Renormalization
		   wind = wind ./ sqrt(qscale);
		   wha = wind .* xhat;      
		   w    = mtlb_ifft(wha);
		   cwt(1:n,kscale) = real(w)';//kscale
		   kscale = kscale+1;
		end
		scale  = scale .*2;
    end 
endfunction