scalogram periods off by factor of sampf

[rickyegeland]

Example code analyzes a simple sine wave with a period of 17.0. scalogram plots power around 170.0

import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
import cwt

# make signal
P = 17.0 # period, s
D = 100   # duration, s
N = 1000 # number of samples
x = np.linspace(0, D, N) # time, s
data = np.sin(2*np.pi/P*x)
label = 'P=%0.1f sine wave' % P

# plot signal
plt.plot(x, data)
plt.title(label)
plt.xlabel('Time [s]')
plt.ylabel('signal')

# continuous wavelet transform
scales = np.arange(200)    
mother = cwt.Morlet(len_signal = len(data), scales = scales)
wavelet = cwt.cwt(data, mother)

# plot with builtin function
wavelet.scalogram()

# plot myself
y = mother.scales / mother.fc / 10.0 # bugfix: divide by 10?
fig = plt.figure()
ax = fig.gca()
contf = ax.contourf(x, y, np.abs(wavelet.coefs)**2)
fig.colorbar(contf, ax=ax, orientation='vertical', format='%2.1f')
ax.set_xlim((x[0], x[-1]))
ax.set_title('scalogram')
ax.set_xlabel('time')
ax.set_ylabel('period')
plt.show()

[rickyegeland]

Actually, "10" is just a result of my sampling frequency. The correction should be:

y = mother.scales / mother.fc / mother.sampf # bugfix

but this of course requires that mother.sampf was properly set, which is not guaranteed. The API is somewhat fragile in that forgetting to set certain values will result in silently producing an incorrect result. An API which calculates all necessary parameters (len_signal, scales, sampf) from an input time series would be better.

[lesserwhirls]

Greetings @rickyegeland!

First off, thank you for the report. It's good to see someone is looking at the code!

Second, when you created the motherwavelet mother, did you set the kwarg sampf? If sampf, the sampling frequency, isn't set, a default of 1 is used.

I agree the API is a bit fragile in terms of forgetting to set values like sampf. It's possible that the mother wavelet could be passed into cwt as a class, and initialized within cwt. That would cut down on one step (creating a mother wavelet), while allowing the signal length to be computed automatically. The scales could be computed automatically based on the COI. However, if the signal is long (for example, the sonic anemometer data I worked with had 36000 - 72000 points), it's nice to be able to specify the particular range of scales that you wish to analyze.

sampf is a bit harder to handle automatically without making some assumptions. If the input time, x, was a numpy array or list of datetime objects, then we could do something smarter there. Of course, maybe I'm thinking about it too much, and we could simply compute sampf from x if it was passed into cwt as well as the mother wavelet.

maybe an api like:

cwt(x, data, mother)

[rickyegeland]

Initially (when I submitted the bug) I wasn't setting sampf. Now I am, however in the current version sampf is not used to properly set the periods in scalogram().

I think the API can be improved, either in the way you suggest, or by implementing cwt() as a method of MotherWavelet. Also, we could allow for the user to set the scales explicitly or allow them to be calculated.

# Ex. 1: explicitly set scales
mother = Morlet()
result = mother.cwt(t, data, scales=myscales)

# Ex. 2: automagically calculate scales
mother = Morlet()
result = mother.cwt(t, data)
# or
result = Morlet().cwt(t, data)
# or
Morlet().cwt(t, data).scalogram()

In any case, parameters like len_signal and sampf should never be necessary as they can be extracted from a time axis, which I expect should always be available to a user. In my current tests I'm setting them up like this.

# t, time axis and data, the signal prepared above...
N = t.size                                                                                                              
T = t[-1] - t[0]                                                                                                        
dt = 1.*T/N                                                                                                             
sampf = 1./dt                                                                                                           
fc = 0.849 # taken from code...                                                                                         
Nscales = int(T * fc * sampf) # compute Nscales to scan from P=0 to P=T                                                                                                                                              
mother = cwt.Morlet(len_signal=N, scales=np.arange(Nscales), sampf=sampf)                                               
wavelet = cwt.cwt(data, mother)

Although as you mention, if a cone of influence is used it doesn't make sense to calculate scales all the way up to T. However, that should be optional too, because if the signal is known to be cyclical there is no need for a cone of influence.

I'm willing to work on these changes if you agree. I would also make a simple unit test and document the usage of this module. I need python wavelet analysis for my research, and although there are alternatives out there, this code is the cleanest implementation I've seen. Maybe if it is improved and tested enough it can make its way into scipy, where I previously wasted many hours trying to figure out how their cwt and Morlet are supposed to work. (They don't. At least not together.)

By the way, I am also working at NCAR (HAO).

[lesserwhirls]

Hi Ricky,

I'd love to work together on this. You are at HAO (Center Green?). Maybe we could get together sometime soon and talk in person about the API.

Thanks!

Sean

On Fri, Jan 30, 2015 at 2:33 PM, Ricky Egeland notifications@github.com wrote:

Initially (when I submitted the bug) I wasn't setting sampf. Now I am, however in the current version sampf is not used to properly set the periods in scalogram().

I think the API can be improved, either in the way you suggest, or by implementing cwt() as a method of MotherWavelet. Also, we could allow for the user to set the scales explicitly or allow them to be calculated.

Ex. 1: explicitly set scales

mother = Morlet() result = mother.cwt(t, data, scales=myscales)

Ex. 2: automagically calculate scales

mother = Morlet() result = mother.cwt(t, data)# or result = Morlet().cwt(t, data)# or Morlet().cwt(t, data).scalogram()

In any case, parameters like len_signal and sampf should never be necessary as they can be extracted from a time axis, which I expect should always be available to a user. In my current tests I'm setting them up like this.

t, time axis and data, the signal prepared above...

N = t.size T = t[-1] - t[0] dt = 1.T/N sampf = 1./dt fc = 0.849 # taken from code... Nscales = int(T * fc \ sampf) # compute Nscales to scan from P=0 to P=T mother = cwt.Morlet(len_signal=N, scales=np.arange(Nscales), sampf=sampf) wavelet = cwt.cwt(data, mother)

Although as you mention, if a cone of influence is used it doesn't make sense to calculate scales all the way up to T. However, that should be optional too, because if the signal is known to be cyclical there is no need for a cone of influence.

I'm willing to work on these changes if you agree. I would also make a simple unit test and document the usage of this module. I need python wavelet analysis for my research, and although there are alternatives out there, this code is the cleanest implementation I've seen. Maybe if it is improved and tested enough it can make its way into scipy, where I previously wasted many hours trying to figure out how their cwt and Morlet are supposed to work. (They don't. At least not together.)

By the way, I am also working at NCAR (HAO).

— Reply to this email directly or view it on GitHub https://github.com/Unidata/pyCWT/issues/1#issuecomment-72274003.