# 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, t[-1], freq[-1], freq], 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, 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.

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

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

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