feature request: discrete wavelet variance estimation

[rgommers]

Originally requested on nigma/pywt#8 by @fairymane :

Hi,
I am trying to implement some wavelet analysis in Python for some R code wrote previously.

One function used in earlier R version code is discrete wavelet variance estimation ( 'wmtsa', page69,http://cran.r-project.org/web/packages/wmtsa/wmtsa.pdf). But I couldn't find equivalent function to implement in pywt package. Does any one knows alternative solution to that if I want to implement the similar function in Python?

Thanks!
Tao

[rgommers]

And my comment:

The DWT based version might not be too hard to implement, but the R default is based on MODWT which is not present in PyWavelets. So that likely won't materialize any time soon. @fairymane do you have an idea about the performance difference between DWT and MODWT, and would DWT-based variance estimation be of use to you?

[grlee77]

I think the new norm argument to swt makes this doable for the undecimated case now. See the new demo that was added.

[jier]

" Hi,
I am trying to implement some wavelet analysis in Python for some R code wrote previously.

One function used in earlier R version code is discrete wavelet variance estimation ( 'wmtsa', page69,http://cran.r-project.org/web/packages/wmtsa/wmtsa.pdf). But I couldn't find equivalent function to implement in pywt package. Does any one knows alternative solution to that if I want to implement the similar function in Python?

Thanks!
Tao"

Any update on this ?? Especially wavelet variance method for the DWT case.

[grlee77]

I think as of the 1.1 release of PyWavelets last year, you can do this kind of analysis with the existing swt functions when you set norm=True. Does the following demo look like what you want?:
https://github.com/PyWavelets/pywt/blob/master/demo/swt_variance.py

(There is some discussion of the feature and figures created from the demo in #476)

[grlee77]

@jier, upon re-reading I see maybe you are specifically requesting this for the standard, decimated DWT rather than SWT?

I think for the DWT case, you just need to rescale the wavelet filterbank by 1/sqrt(2) and then can proceed as otherwise in that demo, only replacing pywt.swt with pywt.wavedec. Here is an explicit example of how to rescale the sym4 wavelet.

import numpy as np
import pywt
sym4 = pywt.Wavelet('sym4')
sym4_normalized = pywt.Wavelet(
    'sym4_normalized',
    filter_bank=[np.asarray(f)/np.sqrt(2) for f in sym4.filter_bank]
)
# Then use `wavelet=sym4_normalized` in the calls to `pywt.wavedec`

The link the the wmtsa PDF above was broken when I tried it, but here is a mirror. The following is a quote from the wavVar docstring from that package which seems to be the function in that library corresonding to this. Note that they recommend to use SWT vs. DWT:

"While DWT wavelet coefficients can also be used, the statistical properties are inferior to those of the MODWT wavelet variance. "

I think in practice you will find substantial variation in the subband variance estimates for various shifts of the input signal because the DWT is not shift-invariant like the SWT.

Another point made in the wmtsa docs is that you may want to exclude some coefficients near the boundaries to avoid bias in the estimates.

[jier]

@grlee77 Yes very much thanks, for the second reply! I indeed need the dwt as I’m working with climate data and I am only using the approximation coefficients for my analysis. Which will not be useful for me to use the swt with norm=True, trim_approx=True, because then I only retain the last approximation coefficient. But maybe is my understanding of wavelet analysis not that really strong.

[jier]

By excluding the coefficients near the boundaries you mean, excluding the begin and end values for each coefficient in each level?

[jier]

By excluding the coefficients near the boundaries you mean, excluding the begin and end values for each coefficient in each level?

I tried what you suggested and it worked!