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.