rwt.sci

来自「小波分解源代码」· SCI 代码 · 共 73 行

function rwt = RWT(x,nvoice,wavelet,oct,scale)
// RWT -- Real Wavelet Transform
//  Usage
//    rwt = RWT(x,nvoice,wavelet)
//  Inputs
//    x        signal, dyadic length n=2^J, real-valued
//    nvoice   number of voices/octave
//    wavelet  string 'Gauss', 'DerGauss','Sombrero', 'Morlet'
//  Outputs
//    rwt      matrix n by nscale where
//             nscale = nvoice .* noctave
//
//  Description
//     see sections 4.3.1 and 4.3.3 of Mallat's book  
//
//  Copyright Aldo I Maalouf
  
  [lhs,rhs]=argn();
	if rhs<4,
		oct = 2;
		scale = 4;
	end	
// preparation
	x = ShapeAsRow(x);
	n = mtlb_length(x);
	xhat = mtlb_fft(x);
	//xhat=xhat';
	xi   = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*%pi/n);
	
// root 
	omega0 = 5;
	
//	noctave = floor(log2(n))-2;
//	noctave = floor(log2(n))-1;
	noctave = floor(log2(n))-oct;
	nscale  = nvoice .* noctave;
	
	rwt = zeros(n,nscale);

	kscale  = 1;
//	scale   = 4;
//	scale = 16;

	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);
		   
		   rwt(1:n,kscale) = real(w)';//kscale
		   kscale = kscale+1;
		end
		scale  = scale .*2;
    end 
    endfunction