PyWavelets
MODWT Transform
Hi all,
It seems to me that there's a problem with the pywt.swt function, possibly with documentation only. The docs say
"Returns list of approximation and details coefficients pairs in order similar to wavedec function:: [(cAn, cDn), ..., (cA2, cD2), (cA1, cD1)]"
which implies that this is actually a multi-resolution analysis, which should mean that the sum $\sum\limits_{k=1}^m cDk + cAm = X$ (X denotes the original signal) for each $m \le n$... this does not happen even for $n=1$, i.e. cD1 + cA1 = X. Anyone know what's going on here?
Are the coefficients actually the wavelet and scaling coefficients? FYI the similar R "wavelets" package function (mra) does satisfy those relations.
Have you checked how this behaves in Matlab's swt? I am checking it right now, but I frankly don't understand Matlab's output structure and it does not seem very well documented.
I have not figured out the format of the output coefficients from Matlab's swt, but at least I can see that if I do a 1-level swt of a vector and then add together the resulting approximation and detail output vectors, the sum does not equal the input vector either.
Yes, the current implementation requires convolving both coefficients again with the filter, and adding the two coefficients together and dividing by 2. The result will be circularly shifted by 2**(level-1), which is easy to correct for.
I have work on an swt branch that implements this behaviour, which matches MATLAB. It conflicts with R's behaviour, and I would very much like to know which is correct. I'm not 100% sure how R's implementation works, and haven't had time to dig into their source/references to work it out.
R uses a 'MODWT' transform, whereas MATLAB uses the "algorithme à trous" to perform shift-invariant transforms. The difference is in how the 'stretched' wavelet filters are calculated.
MATLAB uses something like this, with j being the level, and g the approximation filter. The detail filter h is calculated in the same way.
g[j] = numpy.zeros(len(g[0]) * 2 ** (j - 1))
g[j][::2 ** (j - 1)] = g[0]
MODWT uses different filters (let's say g' and h' for approximation and detail respectively). There is a public introduction available here, and derivation in Percival & Walden, 2000, Wavelet Methods for Time Series Analysis, Equation 95e. Note these refer to g and h as defined above.
g'[j] = convolve(g[0], g[1], ... g[j])
h'[j] = convolve(g'[j], h[j])
I will double check this tomorrow but just wanted to share my findings.
Note that the MODWT filters should allow one to directly compute the jth level DWT (with appropriate downsampling). It doesn't actually reduce the number of operations (the filter is 2**j times longer) but removes the need to allocate scratch space for each level of the transform and probably is cache-friendlier since there is only one pass over the input.
There has been more recent discussion about implementing a MODWT (yay!), with a DFT approach being suggested in #200.
My first instinct would be that this could be slower than a naive convolution, as wavelet filters tend to be small and the corresponding signal quite large. Taking the DFT of the filter and signal, interpolating the filter DFT out to the same length as the signal, doing the per-element multiply, and then the inverse DFT might well be slower than just doing len(filter) * len(signal) operations in the naive version. This is pure speculation, so benchmarks are very much needed.
I agree that FFT-based is not necessarily faster. In a particular case I looked at for DWTs in Matlab in the past, I think FFT-based was faster for filter lengths >16, but obviously the exact cutoff is going to depend on the specific implementation.
