Units of "widths" argument to scipy.signal.cwt() function
Asked Answered
P

1

10

I am confused about the widths parameter that gets passed to scipy.signal.cwt() and by extension to scipy.signal.find_peaks_cwt(). A previous and very helpful Stack Overflow question (and pointers therein) explained away most of my confusion. The widths is an array of scales by which to stretch the wavelet before convolution with your data.

The point which confused me still is, what are the units of the elements of widths? Does a width of 1 mean the wavelet is stretched to be one "index" wide, where index is the distance between elements of data? At first I assumed this was the case, but (a) the widths can take non-integer values, and (b) cwt() results can vary depending on the widths.

Here is some code which illustrates my confusion. Why do the last two lines give different results?

#generating an arbitrary signal with overlapping gaussian peaks with various 
npeaks = 6
support = np.arange(0,1.01,0.01)
pkx = np.array([0.2, 0.3, 0.38, 0.55, 0.65]) #peak locations
pkfun = sum(stats.norm.pdf(support, loc=pkx[i], scale=0.03) for i in range(0,npeaks-1))

#finding peaks for two different setting of widths
pkindsOne = sig.find_peaks_cwt(pkfun, widths = np.arange(4,6), wavelet = sig.ricker)
pkindsTwo = sig.find_peaks_cwt(pkfun, widths = np.arange(4,6.4), wavelet = sig.ricker)

#printing to show difference between calls
for ind, el in enumerate(pkindsTwo):
    print el, pkindsOne[ind]
20 20
36 36
38 38
55 55
63 66
66 91
91

The results are close, but the second call finds one spurious peak at element 63 of the input data. Thus I'm not convinced that the units of widths are indices of the data vector. But what else could they be? If not, what are the units of the widths? cwt() and find_peaks_cwt() never know about or see any x-axis units (e.g. the support vector I define in my code), so what am I missing? When, practically speaking, does it ever make sense to use non-integer widths?

Pontificals answered 3/3, 2015 at 1:4 Comment(0)
B
8

I had the same question myself. Looking at the source code, my best guest is that the units are in "number of samples". The key code line within scipy.signal.wavelets.cwt is:

wavelet_data = wavelet(min(10 * width, len(data)), width)

Here, "wavelet" is a function (builder of the mother wavelet) which receives parameters "length_of_wavelet" and "width_of_wavelet" in number of samples. The reason why widths can still be a non-integer value is that (if I'm not wrong) it represents the scaling factor, which can take any real positive number as it is just a factor of the formulation that affects the shape of the wavelet.

Batwing answered 14/4, 2015 at 10:23 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.