Implementation of iswt
Mike Marino
I just wanted to share some code I devised to make a inverse
stationary wavelet transform using the basis averaging method. The
basic idea is to add iswt to swt, just as there already exists idwt
and dwt, etc. I post this to save others from having to do it again,
but I hope it hasn't already been done in the package and I just
missed it! I haven't speed tested and I haven't extensively tested
(e.g. the code assumes a dyadic array), but hopefully it will save
some work for others. Also, maybe a version could make it into the
package? The code is designed so that the following is true:
>>> output = pywt.swt([0,1,2,3,4,5,6,7], 'haar', level=3)
>>> output
[(array([ 9.89949494, 9.89949494, 9.89949494, 9.89949494,
9.89949494,
9.89949494, 9.89949494, 9.89949494]), array([-5.65685425,
-2.82842712, 0. , 2.82842712, 5.65685425,
2.82842712, 0. , -2.82842712])), (array([ 3., 5.,
7., 9., 11., 9., 7., 5.]), array([-2., -2., -2., -2., -2.,
2., 6., 2.])), (array([ 0.70710678, 2.12132034, 3.53553391,
4.94974747, 6.36396103,
7.77817459, 9.19238816, 4.94974747]), array([-0.70710678,
-0.70710678, -0.70710678, -0.70710678, -0.70710678,
-0.70710678, -0.70710678, 4.94974747]))]
>>> iswt(output, 'haar')
array([ 3.33066907e-16, 1.00000000e+00, 2.00000000e+00,
3.00000000e+00, 4.00000000e+00, 5.00000000e+00,
6.00000000e+00, 7.00000000e+00])
Cheers, Mike
***** Code follows:
import math
import pywt
import numpy
def iswt(coefficients, wavelet):
"""
Input parameters:
coefficients
approx and detail coefficients, arranged in level value
exactly as output from swt:
e.g. [(cA1, cD1), (cA2, cD2), ..., (cAn, cDn)]
wavelet
Either the name of a wavelet or a Wavelet object
"""
output = coefficients[0][0].copy() # Avoid modification of input
data
#num_levels, equivalent to the decomposition level, n
num_levels = len(coefficients)
for j in range(num_levels,0,-1):
step_size = int(math.pow(2, j-1))
last_index = step_size
_, cD = coefficients[num_levels - j]
for first in range(last_index): # 0 to last_index - 1
# Getting the indices that we will transform
indices = numpy.arange(first, len(cD), step_size)
# select the even indices
even_indices = indices[0::2]
# select the odd indices
odd_indices = indices[1::2]
# perform the inverse dwt on the selected indices,
# making sure to use periodic boundary conditions
x1 = pywt.idwt(output[even_indices], cD[even_indices],
wavelet, 'per')
x2 = pywt.idwt(output[odd_indices], cD[odd_indices],
wavelet, 'per')
# perform a circular shift right
x2 = numpy.roll(x2, 1)
# average and insert into the correct indices
output[indices] = (x1 + x2)/2.
return output
Nicolas Barbey
Thanks for your contribution. I am interested into this too. I briefly
tested this code and it seems to work for me in quite some cases.
But what I really need is the 2D counterpart so I consider writing an
implementation of the "iswt2" function. Any advice on this one ?
Also, is there any chance that it would make into PyWavelets trunk ?
On this line of though, is there any coding practice or style that you
would require for a contribution to this package ?
As I am fairly new to wavelets it may take some time for me to code
this but I am willing to do my best :)
Cheers,
Nicolas
indi...@gmail.com
Sorry for replying to a 3 year old email thread.
I was looking for iswt in pywt, but couldn't find it. Is it implemented in pywt?
If not, i think i shall go ahead with Marino's algorithm.
Thank you.
-cheers
joe
Mike Marino
--
You received this message because you are subscribed to the Google Groups "PyWavelets" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pywavelets+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
Mike Marino
Sure, no problem.
Cheers,
Mike
Hi Mike,
Better late than never, I hope, we are currently working on including iswt and iswt2 in PyWavelets (https://github.com/PyWavelets/pywt/issues/57#issuecomment-126419900). I was hoping to use your proposed iswt here for that. Would that be OK by you?
Best,
Thomas