[SciPy-dev] Perhaps an addition to scipy.signal.wavelets?

[iCy-fLaME]

I been messing around with the continuous wavelet transform lately,
but not seems to be able to find much resources in scipy about
wavelets.

The Morlet wavelet is a pretty standard one to use, so I coded this
function to generate the wavelet, hope it can be some use for others.
Might require some modification, because I coded this as a standalone
function.


########
# Morlet #
########

def morlet(M, w = 5.0, s = 1.0, cplt = True):
	"""
	Returns the complex Morlet wavelet with length of M
	
	Inputs:
		M	Length of the wavelet
		w	Omega0
		s	Scaling factor, windowed from -s*2*pi to +s*2*pi
		cplt	Use the complete or standard version
	
	The standard version:
	pi**-0.25 * exp(1j*w*x) * exp(-0.5*(x**2))
	
	Often referred to as simply the Morlet wavelet in many text,
		also commonly use in practice. However, this simplified version
		can cause admissibility problem at low w.
	
	The complete version:
	pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2))
	
	Complete version of the Morlet wavelet, with the correction
		term to improve admissibility. For w is greater than 5, the
		correction term will be negligible.
	
	The energy of the return wavelet is NOT normalised according to s.
	
	NB: The fundamental frequency of this wavelet in Hz is given by:
		
		f = 2*s*w*r / M
		
		r - Sampling Rate
	
	"""
	
	from scipy import linspace
	from scipy import pi
	from scipy import exp
	from scipy import zeros
	
	correction = exp(-0.5*(w**2))
	c1 = 1j*w
	s *= 2*pi
	
	xlist = linspace(-s, s, M)
	output = zeros(M).astype('F')
	
	if cplt == True:
		for i in range(M):
			x = xlist[i]
			output[i] = (exp(c1*x) - correction) * exp(-0.5*(x**2))
	else:
		for i in range(M):
			x = xlist[i]
			output[i] = exp(c1*x) * exp(-0.5*(x**2))
	
	output *= pi**-0.25
	
	return output

#################################################


-- 
iCy-fLaME
The body maybe wounded, but it is the mind that hurts.