# Continuous wavelet transform (CWT)

Hello! Please tell me if there is a way to build a CWT (or to be more precise, a CWT analog for a discrete signal) of a signal in Sage. The task is to obtain a time evolution of the spectrum (more or less) for the recorded signal in the form of time series. I've found the good examples of what I mean here.

Thanks.

edit retag close merge delete

Sort by » oldest newest most voted

A self-answer: I downloaded an mlpy (Machine Learning PYthon) package and installed it (I use Gentoo, so the system-widely installed python modules are visible in Sage). After that I applied its functions to my time series in the following way:

signal01_02_array = numpy.array(signal01_02_AC) # 'signal01_02_AC' is a regular 1D python list containing my time series

import matplotlib.pyplot as plt
import mlpy

omega0 = 8
spec, scale = mlpy.cwt(signal01_02_array, dt=0.001, dj=0.05, wf='morlet', p=omega0,
extmethod='zeros', extlength='powerof2') # 'dt' is a time step in the time series

freq = scale[1:] # building a 'Frequency' scale (this is absolutely NOT correct, because the scale <-> frequency correspondence strongly depends on the certain wavelet form)

t = numpy.array(signal01_02_array) #
# building a 'Time' scale
for i in range(0,len(t)):          #
t[i] = float(i)/1000           #

fig = plt.figure() # creating an empty Matplotlib figure
ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2], xticks=[0,1,2,3,4,5,6,7,8,9]) # axis for the initial signal
ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], xticks=[0,1,2,3,4,5,6,7,8,9], xlabel="Time, s", ylabel="Scales") # axis for the resulting wavelet transform data
ax3 = fig.add_axes([0.83,0.1,0.03,0.6], xlabel="Magnitude") # axis for a colour bar
ax1.plot(t, signal01_02_array, 'k') # plotting the initial time series

img = ax2.imshow(numpy.abs(spec), extent=[t[0], t[-1], freq[-1], freq[0]], aspect='auto') # plotting a CWT data
fig.colorbar(img, cax=ax3) # building a colour bar basing on the colours present in CWT data image
plt.show()


This works perfectly fine in a pure Python-2.7 shell, but fails in Sage with the following error:

Traceback (most recent call last):                           extmethod='none', extlength='powerof2')
File "", line 1, in <module>

File "/tmp/tmpVOEsSJ/___code___.py", line 9, in <module>
extmethod='none', extlength='powerof2')
File "/usr/lib64/python2.7/site-packages/mlpy/_cwt.py", line 160, in cwt
s = scales(x.shape[0], dj, dt, s0)
File "/usr/lib64/python2.7/site-packages/mlpy/_cwt.py", line 76, in scales
s[i] = s0 * 2**(i * dj)
TypeError: unsupported operand type(s) for ** or pow(): 'int' and 'sage.rings.real_mpfr.RealNumber'


When 'dj' is a whole number (like 1, 2, ...), everything works fine in Sage too.

more

There is a way to overcome the problem mentioned above. One needs to typecast 'dj' to float. Like this:

spec, scale = mlpy.cwt(signal01_02_array, dt=0.001, dj=float(0.05), wf='morlet', p=omega0,
extmethod='zeros', extlength='powerof2') # 'dt' is a time step in the time series


and it will work in Sage also.

more

I think I can address your last problem: 2 is being treated like a python int, but the exponent, apparently, is a sage.rings.real_mpfr.RealNumber. I can reproduce this error with

sage: int(2)**.5
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
...
TypeError: unsupported operand type(s) for ** or pow(): 'int' and 'sage.rings.real_mpfr.RealLiteral'


On the other hand, both of the following work:

sage: x = Integer(2)**.5; x
1.41421356237310
sage: type(x)
<type 'sage.rings.real_mpfr.RealNumber'>

sage: y = int(2)**float(.5); y
1.4142135623730951
sage: type(y)
<type 'float'>


So depending on what you need, you could either change dj to a python float, or change 2 to a sage integer.

more

Yes, I think so too, but the integer '2' is inside one of the mlpy files. Probably it is necessary to patch this file during installation to make it compatible with Sage.

( 2011-06-22 17:21:33 +0100 )edit

Ah; I wondered if that was the case. Can you simply make dj a float to work around the problem?

( 2011-06-22 17:33:22 +0100 )edit

Well... surprisingly (for me) it works! :) I'll post one more answer to explicitly illustrate this solution. Thanks!

( 2011-06-22 17:47:29 +0100 )edit

Thought this looked familiar. Fixing int-to-RealNumber is http://trac.sagemath.org/sage_trac/ticket/10736, and it has a positive review.

( 2011-06-22 20:33:50 +0100 )edit