PyWavelets / pywt

PyWavelets - Wavelet Transforms in Python
http://pywavelets.readthedocs.org
MIT License
2.07k stars 482 forks source link

CDF 9/7 wavelet #176

Closed nareto closed 8 years ago

nareto commented 8 years ago

I would like to use the CDF 9/7 wavelets on a 1D signal. This is the test code I have, which doesn't work (i.e. the reconstructed signal is different, though somewhat similar....)

import numpy as np
import pywt

factor = 1
cdf97_an_lo = factor*np.array([0.026748757411, -0.016864118443, -0.078223266529, 0.266864118443,\
                        0.602949018236, 0.266864118443, -0.078223266529,-0.016864118443, \
                        0.026748757411])
cdf97_an_hi = factor*np.array([0, 0.091271763114, -0.057543526229,-0.591271763114,1.11508705,\
                        -0.591271763114,-0.057543526229,0.091271763114,0 ])
cdf97_syn_lo = factor*np.array([0,-0.091271763114,-0.057543526229,0.591271763114,1.11508705,\
                         0.591271763114 ,-0.057543526229,-0.091271763114,0])
cdf97_syn_hi = factor*np.array([0.026748757411,0.016864118443,-0.078223266529,-0.266864118443,\
                         0.602949018236,-0.266864118443,-0.078223266529,0.016864118443,\
                         0.026748757411])
cdf97 = pywt.Wavelet('cdf97', [cdf97_an_lo,cdf97_an_hi,cdf97_syn_lo,cdf97_syn_hi])
wav = cdf97
mode = 'periodization'
sig= np.arange(8)
a,d = pywt.dwt(sig,wav,mode)
print("approx: %s \n detail: %s \n" %(a,d))
rec_sig = pywt.idwt(a,d,wav,mode)
print(rec_sig)

the output is:

  approx: [ 0.45329098  3.21399006  4.92092289  5.41179607] 
  detail: [  2.69825893e-01  -7.37399958e-09  -7.30174117e-01   4.46034819e+00] 

 [-0.2161735   1.2161735   3.54925496  4.54925498  4.21617348  4.78382651
   8.45074501  1.45074503]

am I doing something wrong? I'm kind of new to wavelet theory. I noticed however that once I create the pywt.Wavelet object, a 0 gets added to all the filters, making them of length 10 (while they should be of length 9) - is this expected behaviour? I.e.

print(cdf97.get_filters_coeffs()[0])

prints:

 [0.026748757411,
  -0.016864118443,
  -0.078223266529,
  0.266864118443,
  0.602949018236,
  0.266864118443,
  -0.078223266529,
  -0.016864118443,
  0.026748757411,
  0.0]
grlee77 commented 8 years ago

I am not sure about the padding without looking into it more closely, but the built-in bior4.4 wavelet is the same as CDF 9/7 aside from a different normalization of the coefficients. (e.g. dec_lo is scaled by sqrt(2) and dec_hi by 1/sqrt(2)).

Based on comparing where the 0's are in the bior4.4 coefficients, the following slight modification to your cdf97 coefficients to have the same symmetry gives the expected exact reconstruction (to within numerical precision). Note that two of the filters have a zero appended and two have it prepended.

factor=1
cdf97_an_lo = factor*np.array([0, 0.026748757411, -0.016864118443, -0.078223266529, 0.266864118443,\
                        0.602949018236, 0.266864118443, -0.078223266529,-0.016864118443, \
                        0.026748757411])
cdf97_an_hi = factor*np.array([0, 0.091271763114, -0.057543526229,-0.591271763114,1.11508705,\
                        -0.591271763114,-0.057543526229,0.091271763114,0, 0 ])
cdf97_syn_lo = factor*np.array([0, -0.091271763114,-0.057543526229,0.591271763114,1.11508705,\
                         0.591271763114 ,-0.057543526229,-0.091271763114,0, 0])
cdf97_syn_hi = factor*np.array([0, 0.026748757411,0.016864118443,-0.078223266529,-0.266864118443,\
                         0.602949018236,-0.266864118443,-0.078223266529,0.016864118443,\
                         0.026748757411])
cdf97 = pywt.Wavelet('cdf97', [cdf97_an_lo,cdf97_an_hi,cdf97_syn_lo,cdf97_syn_hi])
nareto commented 8 years ago

thank you, it seemed strange these filters weren't present in pywt. Now I'm using bior4.4 and reconstruction works as expected