regeirk/pycwt

A Python module for continuous wavelet spectral analysis. It includes a collection of routines for wavelet transform and statistical analysis via FFT algorithm. In addition, the module also includes cross-wavelet transforms, wavelet coherence tests and sample scripts.

pycwt's Issues

Deprecated function in scipy

When using wavelet.cwt I get the following warning about a future deprecation in scipy:

/home/jorge/.local/lib/python3.5/site-packages/scipy/fftpack/basic.py:159: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use arr[tuple(seq)] instead of arr[seq]. In the future this will be interpreted as an array index, arr[np.array(seq)], which will result either in an error or a different result.
z[index] = x

calculation of sig2 in sample_xwt.py has a mistake

Same data1 is used for sig2 instead of data2

current code
sig2 = wavelet.significance(1.0, dt, cwt2[1], 0, data1['alpha'], wavelet=mother)
should be:
sig2 = wavelet.significance(1.0, dt, cwt2[1], 0, data2['alpha'], wavelet=mother)

BSD Licence

Hey,

I really want to get this wavelet code licensed so that maybe it can be added to SciPy. They require a BSD license for their code if I recall.

I was wondering what code you translated into python and what you wrote yourself?
If I can contact everyone that you translated code from and get their permission or for them to licence their own code under BSD. We should be able to licence this module with BSD.

Correlations package not found

Hi,

I do not manage to import pycwt since the correlations package for the acvf function is not found. I cannot find this package online. Can sombody help me on this?

Cheers,
Heiner

Python3 compatability

To get pycwt to work on python 3.5.1, I had to adjust wavelet.py, by commenting out these lines
if isinstance(wavelet, unicode) or isinstance(wavelet, str):
wavelet = mothers[wavelet]()

I understand they are supposed to call the constructor of the correct pycwt.mother.wavelet class, but the dictionary lookup seems to be breaking over a problem to do with unicode.

Inconsistent results?

I copy the code to notebook on my machine and tested.
The result as this png:

As the official result is:

pycwt document

The biggest difference is the minimum value of the original data from document is letter than -2, but my test is greater than -2? Is that why they had inconsistent results?

IndexError: too many indices -- in wavelet.wct code

Hello,

When I run your sample_xwt.py code, I get an error when it calculates the wavelet coherence (WTC). I see the error running the latest Enthought Canopy on both OSX and Windows.

It seems to complete the coherence significance calculation, it reaches 100%, however is followed by the error:

Calculating wavelet coherence significance... 100.0%, 6m03.4s (0.0s, 1.1s) 
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "wav.py", line 953, in wct
    sig = wct_significance(a1, a2, significance_level=0.95, **kwargs)
  File "wav.py", line 1072, in wct_significance
    sig95[s] = interp(significance_level, P, R2y[s:, sel])
IndexError: too many indices

Did you modify something to get it to run on your end?

inverse cwt

I believe there is a square root missing (sj) in the calculation of the inverse wavelet transform. It should be:

iW = (dj * np.sqrt(dt) / wavelet.cdelta * wavelet.psi(0) *
(np.real(W) / np.sqrt(sj)).sum(axis=0))

Preferred citation

Hi, thank you for making and maintaining this library! What is the preferred way to cite pycwt?

Biorthogonal Wavelets

Hello!

I've used your package and it works great. I use it for ECG Denoising an Delineation. In articles researchers often use Biorthogonal wavelets, such as 'bior1.5' (for example in Vítek M. et all (2010) - http://bs2010.biosignal.cz/papers/1042.pdf). I'd like to add such mother wavelet to your code, but unfortunately I can't found any formulae for it. Could you help me with it, please?

typo in wavelet.py

top line, looks like an "int" accidentally got pasted in PR #43

Can't not import pyssim

I want to use this algorithm in my project, but I found that this module cannot be imported using 'from pyssim import ssim'. Can someone tell me how to resolve this?

red noise generator may not function as intended

I think the rednoise function does not produce the results desired. Here is a minimal example based on the function.

from scipy.signal import lfilter
import matplotlib.pyplot as plt
import numpy as np


# This should be strongly non-white when g is close to 1
g = 0.9
# standard deviation of the gaussian white noise
a = 1
# Just a short sample to make visual assessment simpler 
N = 100
tau = int(np.ceil(-2 / np.log(np.abs(g))))

# Save the noise for comparison
noise = np.random.randn(N + tau, 1) * a
# Apply the filter
yr = lfilter([1, 0], [1, -g], noise)
yr = yr[tau:]

# Compare the two
plt.plot(yr)
plt.plot(noise[tau:])

which produces the attached plot. For some reason the parameters passed to lfilter do not modify the noisy time series in the way that I would have anticipated. I verified that rednoise returns flat spectra, consistent with the white noise being unaltered as in the example.

Am I implementing the rednoise function incorrectly?

AR1 Warning: Cannot place an upperbound on the unbiased AR(1). Series is too short or trend is to large.

When I want to replace the nan values of the dataset with mean value, sometimes it will raise the warning. I can't tell the difference between the problematic data and the normal data. If I use zero or mean/2, this warning will not raise. Does anyone meet this issue?

Release on pip

Hi,
I am using your code to perform CWT on some data of mine, and I need to install it on a bunch of machines. Getting the code from github is not very practical...

Can you release pycwt as a python package in pypi? That way it would be much easier to install, with:

> pip install pycwt

and it can also be included as a dependency by other packages.

Here is a nice howto:

http://peterdowns.com/posts/first-time-with-pypi.html

Thanks!

Rectification scheme

Hello everybody. Does anyone know how can we implement the Liu et al. (2007) rectification scheme in the sample.py code? It is only passed in the plot.py module and I cannot call it in sample.py

ImportError: No module named wavelet

Hello,
I'm trying to run your sample code and am getting stuck when you try to load the wavelet module. I am using enthought canopy.
I can get to "wavelets" (with an 's') through: from scipy.signals import wavelets
However "wavelet" doesn't seem to exist.
Is it your own code? Any help is appreciated!
Cheers,

should not raise Warning ,but use the warning module in python

Raise interrupts the execution while the original matlab ar1 method did not have that intention.

In wav.py in ar1 function

   ``` raise Warning ('Cannot place an upperbound on the unbiased AR(1). '
        'Series is too short or trend is to large.')```

should be replace as below with the warning module function warn()

if D > 0:
        g = (-B - D**0.5) / (2 * A)
    else:
        warnings.warn('Cannot place an upperbound on the unbiased AR(1). '
            'Series is too short or trend is to large.', UserWarning)

Mistake in the red-noise spectrum

pycwt/pycwt/wavelet.py

Lines 252 to 257 in 1bef44e

	# Theoretical discrete Fourier power spectrum of the noise signal
	# following Gilman et al. (1963) and Torrence and Compo (1998),
	# equation 16.
	def pk(k, a, N):
	return (1 - a ** 2) / (1 + a ** 2 - 2 * a * np.cos(2 * np.pi * k / N))
	fft_theor = pk(freq, alpha, n0)

According to equation 16 from Torrence and Compo (1998) this code should look like this

def pk(k, a, N): 
    return (1 - a ** 2) / (1 + a ** 2 - 2 * a * np.cos(2 * np.pi * k / N)) 
fft_theor = pk(freq, alpha, 1) 

This bug leads to different results of usage of significance function with signal and its variance.

Pip install broken

Hi,

the pip install seems to fail (for me) seems to fail with the read("readme.md") -> file not found exception.

great package!

Error: 'Cannot place an upperbound on the unbiased AR(1). '

Hello!
I'm trying to execute sample_XWT.py (inside sample foder) with my own data and all the time this Warning is raising. It happened because, for the second time series (s2), D-variable (within 'ar1' function) is negative.
Does anyone know what couId try?
s2 is a nx1 numpy array (n=58) with positive float numbers.
Thank you very much in advance!

About the time granularity.

Hello, thank you very much for your contribution to enable me to use Python for wavelet analysis. I have a question that I would like to ask, and I wonder if you are comfortable answering it.

In your example data, the time granularity is seasonal, i.e. 0.25 years, and it is continuous full year data. Now I have a time series with a daily granularity, but only observations are available from 1 February to 30 November each year. There are two ways I can think of to do wavelet analysis on such data.

1. fill all the days of the year for which there is no data with zeros, and then perform a wavelet analysis.
2. Define the number of days in the year and define the period for which data is actually available as a year, i.e. 303 days. Afterwards perform a wavelet analysis and apply the definition of a year to all cycles of the analysis.

My personal preference is to use the second approach, but is this theoretically feasible?

raise a error.

while using the function cwt,
raise this line:C:\Anaconda2\lib\site-packages\pycwt-0.2-py2.7.egg\pycwt\wavelet.py:121: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

Does it make some mistake?

Long term goals

I want to start a discussion about what the goals should be to "finish" this module and get it ready for a real release.

To that end, it would be nice to have a list of current problems or improvements/features that users would like to see fixed or added.

These are what I have in mind.

• Tests.
• Expanded mother selection.
• Speeding up of cross-wavelet significance code.
• Taking account of power spectrum bias.
• Adding emails regarding the licence change.
• Documentation.

Anyone else have any suggestions?

Reconstructed waveform from icwt is still different from the original timeseries.

After fixing the bug from PR( #29 ), I test the following code and plot the original timeseries and iwave, but they are not fully same.

W, scales, freqs, coi, _, _ = wavelet.cwt(s, dt, dj, s0, J, mother)
iwave = wavelet.icwt(W, scales, dt, dj, mother)

And it seems in the example of https://pycwt.readthedocs.io/en/latest/, there is the same problem in fig a).
Have anyone met this and how to solve it?
Thanks a lot!

the dof for the global wavelet spectrum will change if scales are in days instead of years?

According to the tutorial sample.py, for the global wavelet spectrum the input parameter "dof" for "wavelet.significance(...,dof=dof)" (lines 91-92) is calculated using dof = N - scales, where N is the size of the input data and the parameter scales has units of years, then this "dof" value is used as "na" in eq (23) of Torrence and Compo (1998). However, if I want to use scales in days for the same problem the dof will change and will be negative when scales > N (N=504 for Nino3 example), which is not consistent.

According to Torrence and Compo (1998) (page 73), they mention that for the global wavelet if the points going into the average of eq (22) are within the cone of influence, then the number of effective points "na" of eq (23) is reduced
by approximately one-half of the number within the COI to reflect the decreased amplitude (and information) within that region. So it means that the effective number of points "na" of eq (23) should be: N - the half of the points inside the COI at each scale? or it might be better to use N - all the points inside the COI at each scale? If so, when we are at the limit of the larger scales, in which all the points are almost inside the COI, na=N-N/2=N/2 for the former case and na=N-N=0 and for the latter one.

In conclusion, what approximation should we use to get the most realistic "na" for eq (23) that does not depend on the units of the parameter scales?

Thanks.
Jeancarlo Fajardo U.

wtc Morlet "boxcar" smoothing on non-linear (2^j) spaced scales

First of all, thank you for this package. I'm glad this repository has come back to life.
I'd like to use this rise in activity as an opportunity to ask a question and/or possibly point to a problem in the implementation.

During wtc calculation (using Morlet wavelet) the wavelet transforms are smoothed in time and scale as suggested in Torrence and Webster (1999). However, this "boxcar" filter in pycwt in scale is done on 2^j spaced scales, which is hardly a boxcar filter. Torrence and Compo (1998) hints at doing something like that for general smoothing in scales over some band, but an average of powers weighted over scales is used. So how should I interpret this implementation of this scale smoothing operation?

I've also asked this question some time ago on DSP SE.

sample.py does not reproduce example graphic

I made a bunch of fixes (see pull #4) so both samples run. The graphics produced by sample_xwt.py appear identical to the graphics provided as example output. However, sampl.py produces a graphic similar to the sample graphic but not identical. The difference should probably be investigated before one trusts the code.

Conflict with matplotlib 1.4.3

I tried to run sample.py and I got the following error:

.../site-packages/matplotlib/init.py:855: UserWarning: text.fontsize is deprecated and replaced with font.size; please use the latter.

I assume there is a conflict with the newer versions of matplotlib. Can there be an easy fix?

Here is a similar issue: mrtommyb/K2fov#3

Rectification of the Bias in the Wavelet Power Spectrum?

I see the power spectrum correction mentioned in the documentation, but when I add the following code, the result is not consistent with the original result, can you please see what is going on?

with rectification:

without rectification :

Below is the original paper pitcure(upper is without rectification):

icwt can't produce the same waveform

import pycwt as wavelet
dt=0.002
mother = wavelet.DOG()
coefs, scales, freqs, coi, fft, fftfreqs = wavelet.cwt(wave, dt,
                                                      mother,freqs=np.arange(2,10,0.2))
iwave = wavelet.icwt(coefs, scales, dt, scales[0]-scales[1], mother)

The waveform changed a lot

Wrong default parameter in function 'wct_significance'

I think there is an error for the parameter 'wavelet' in the function 'wct_significance'. The docstring says it should be an instance of a class, but a string is given.

Function's parameter:

pycwt/pycwt/wavelet.py

Line 532 in d53dbe4

wavelet='morlet', mc_count=300, progress=True,

Docstring:

pycwt/pycwt/wavelet.py

Line 557 in d53dbe4

wavelet : instance of a wavelet class, optional

Edit: I should add that the function '_check_parameter_wavelet' is not used in 'wct_significance'. If one were to use this function with the default parameter, they would get an error. In fact, in the function 'cwt', the function '_check_parameter_wavelet' is used before using 'wct_significance'.

pycwt/pycwt/wavelet.py

Line 472 in d53dbe4

wavelet = _check_parameter_wavelet(wavelet)

pycwt/pycwt/wavelet.py

Line 522 in d53dbe4

sig = wct_significance(a1, a2, dt=dt, dj=dj, s0=s0, J=J,

Monsoon.dat file

I have seen the monsoon.dat file in pycwt. Can anyone explain the file why there are negative values in the data? I have tried to subtract the mean from the values but still I am unable to figure it out.

What should I put in "sample spacing" `dt` parameter for `wct()`?

I am trying to recreate some stuff from a paper, and one point it requires wavelet coherence from matlab to find the symmetry of face. matlab wcoherence() doesn't have this "sample spacing" parameter but pycwt has dt for sample spacing for wct(). I am sorry I am new to signal processing