Ask Your Question
0

Continuous wavelet transform (CWT)

asked 2011-06-21 17:57:03 +0100

v_2e gravatar image

updated 2015-01-14 14:27:16 +0100

FrédéricC gravatar image

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 flag offensive close merge delete

3 Answers

Sort by » oldest newest most voted
1

answered 2011-06-22 16:46:47 +0100

v_2e gravatar image

updated 2011-06-23 01:38:43 +0100

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.

edit flag offensive delete link more
0

answered 2011-06-22 17:51:41 +0100

v_2e gravatar image

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.

edit flag offensive delete link more
0

answered 2011-06-22 17:17:47 +0100

niles gravatar image

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.

edit flag offensive delete link more

Comments

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.

v_2e gravatar imagev_2e ( 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?

niles gravatar imageniles ( 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!

v_2e gravatar imagev_2e ( 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.

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

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-06-21 17:57:03 +0100

Seen: 1,782 times

Last updated: Jun 23 '11