Source code :: cwt
[Return]
[Download]## This code is written by Davide Albanese, <albanese@fbk.eu> and
## Marco Chierici, <chierici@fbk.eu>.
## (C) 2008 Fondazione Bruno Kessler - Via Santa Croce 77, 38100 Trento, ITALY.
## See: Practical Guide to Wavelet Analysis - C. Torrence and G. P. Compo.
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
from numpy import *
import _extend
__all__ = ["cwt", "icwt", "angularfreq", "scales", "compute_s0"]
def angularfreq(N, dt):
"""Compute angular frequencies.
:Parameters:
N : integer
number of data samples
dt : float
time step
:Returns:
angular frequencies : 1d numpy array
"""
# See (5) at page 64.
N2 = N / 2.0
w = empty(N)
for i in range(w.shape[0]):
if i <= N2:
w[i] = (2 * pi * i) / (N * dt)
else:
w[i] = (2 * pi * (i - N)) / (N * dt)
return w
def scales(N, dj, dt, s0):
"""Compute scales.
:Parameters:
N : integer
number of data samples
dj : float
scale resolution
dt : float
time step
:Returns:
scales : 1d numpy array
scales
"""
# See (9) and (10) at page 67.
J = floor(dj**-1 * log2((N * dt) / s0))
s = empty(J + 1)
for i in range(s.shape[0]):
s[i] = s0 * 2**(i * dj)
return s
def compute_s0(dt, p, wf):
"""Compute s0.
:Parameters:
dt :float
time step
p : float
omega0 ('morlet') or order ('paul', 'dog')
wf : string
wavelet function ('morlet', 'paul', 'dog')
:Returns:
s0 : float
"""
if wf == "dog":
return (dt * sqrt(p + 0.5)) / pi
elif wf == "paul":
return (dt * ((2 * p) + 1)) / (2 * pi)
elif wf == "morlet":
return (dt * (p + sqrt(2 + p**2))) / (2 * pi)
else:
raise ValueError("wavelet '%s' is not available" % wf)
def cwt(x, dt, dj, wf="dog", p=2, extmethod='none', extlength='powerof2',res = 2000, fmin = 0, fmax = inf):
"""Continuous Wavelet Tranform.
:Parameters:
x : 1d numpy array
data
dt : float
time step
dj : float
scale resolution (smaller values of dj give finer resolution)
wf : string ('morlet', 'paul', 'dog')
wavelet function
p : float
wavelet function parameter
extmethod : string ('none', 'reflection', 'periodic', 'zeros')
indicates which extension method to use
extlength : string ('powerof2', 'double')
indicates how to determinate the length of the extended data
:Returns:
(X, scales) : (2d numpy array complex, 1d numpy array float)
transformed data, scales
Example:
>>> import numpy as np
>>> import mlpy
>>> x = np.array([1,2,3,4,3,2,1,0])
>>> mlpy.cwt(x=x, dt=1, dj=2, wf='dog', p=2)
(array([[ -4.66713159e-02 -6.66133815e-16j,
-3.05311332e-16 +2.77555756e-16j,
4.66713159e-02 +1.38777878e-16j,
6.94959463e-01 -8.60422844e-16j,
4.66713159e-02 +6.66133815e-16j,
3.05311332e-16 -2.77555756e-16j,
-4.66713159e-02 -1.38777878e-16j,
-6.94959463e-01 +8.60422844e-16j],
[ -2.66685280e+00 +2.44249065e-15j,
-1.77635684e-15 -4.44089210e-16j,
2.66685280e+00 -3.10862447e-15j,
3.77202823e+00 -8.88178420e-16j,
2.66685280e+00 -2.44249065e-15j,
1.77635684e-15 +4.44089210e-16j,
-2.66685280e+00 +3.10862447e-15j,
-3.77202823e+00 +8.88178420e-16j]]), array([ 0.50329212, 2.01316848]))
"""
#x -= mean(x)
lenght = x.shape[0]
if extmethod != 'none':
x = _extend.extend(x, method=extmethod, length=extlength)
w = angularfreq(x.shape[0], dt)
s0 = compute_s0(dt, p, wf)
s = scales(lenght, dj, dt, s0)
freq = (p + sqrt(2.0 + p**2))/(4*pi * s)
ind = where((freq>fmin)&(freq<fmax))
s = s[ind]
x = fft.rfft(x, axis=0)
step = max(int(lenght/res),1)
spec = zeros((len(s),len(w[0:lenght:step]) ),dtype=complex)
stmp = zeros(1)
#wavelet = zeros((len(w)),dtype=complex )
wft = zeros((len(w)),dtype=complex )
for i in range(len(s)):
interv = where((abs(s[i]*w[0:len(w)/2]-p))<3)
wavelet = (1+sign(w[interv]))*exp(-(s[i]*w[interv]-p)**2/2)*sqrt(abs(s[i])/dt)
#wavelet[interv] = waveletb.morletft(stmp, w[interv], p, dt, norm = True)
wft[interv] = x[interv]*wavelet
wft = fft.ifft(wft )
spec[i,:] = wft[0:lenght:step]
wft[:] = 0
return spec, s
def icwt(X, dt, dj, wf = "dog", p = 2, recf = True):
"""Inverse Continuous Wavelet Tranform.
:Parameters:
X : 1d numpy array
transformed data
dt : float
time step
dj : float
scale resolution (smaller values of dj give finer resolution)
wf : string ('morlet', 'paul', 'dog')
wavelet function
p : float
wavelet function parameter
* morlet : 2, 4, 6
* paul : 2, 4, 6
* dog : 2, 6, 10
recf : bool
use the reconstruction factor (:math:`C_{\delta} \Psi_0(0)`)
:Returns:
x : 1d numpy array
Example:
>>> import numpy as np
>>> import mlpy
>>> X = np.array([[ -4.66713159e-02 -6.66133815e-16j,
... -3.05311332e-16 +2.77555756e-16j,
... 4.66713159e-02 +1.38777878e-16j,
... 6.94959463e-01 -8.60422844e-16j,
... 4.66713159e-02 +6.66133815e-16j,
... 3.05311332e-16 -2.77555756e-16j,
... -4.66713159e-02 -1.38777878e-16j,
... -6.94959463e-01 +8.60422844e-16j],
... [ -2.66685280e+00 +2.44249065e-15j,
... -1.77635684e-15 -4.44089210e-16j,
... 2.66685280e+00 -3.10862447e-15j,
... 3.77202823e+00 -8.88178420e-16j,
... 2.66685280e+00 -2.44249065e-15j,
... 1.77635684e-15 +4.44089210e-16j,
... -2.66685280e+00 +3.10862447e-15j,
... -3.77202823e+00 +8.88178420e-16j]])
>>> mlpy.icwt(X=X, dt=1, dj=2, wf='dog', p=2)
array([ -1.24078928e+00, -1.07301771e-15, 1.24078928e+00,
2.32044753e+00, 1.24078928e+00, 1.07301771e-15,
-1.24078928e+00, -2.32044753e+00])
"""
rf = 1.0
if recf == True:
if wf == "dog" and p == 2:
rf = 3.13568
if wf == "dog" and p == 6:
rf = 1.70508
if wf == "dog" and p == 10:
rf = 1.30445
if wf == "paul" and p == 2:
rf = 2.08652
if wf == "paul" and p == 4:
rf = 1.22253
if wf == "paul" and p == 6:
rf = 0.89730
if wf == "morlet" and p == 2:
rf = 2.54558
if wf == "morlet" and p == 4:
rf = 0.92079
if wf == "morlet" and p == 6:
rf = 0.58470
s0 = compute_s0(dt, p, wf)
s = scales(X.shape[1], dj, dt, s0)
# See (11), (13) at page 68
XCOPY = empty_like(X)
for i in range(s.shape[0]):
XCOPY[i] = X[i] / sqrt(s[i])
x = dj * dt **0.5 * sum(real(XCOPY), axis = 0) / rf
return x[Return]