python - CDF 9/7 Discrete Wavelet Transform (Convolution) -

i'm trying write simple self-contained program single level of discrete wavelet transform on 1d list, using cdf 9/7 wavelets, , reconstructs it. i'm using convolution/filter-bank method grasp on how works. in other words, convolve list filter scale coefficients, convolve list different filter wavelet coefficients, starting every other element. upsample (i.e. add zeroes in between elements), apply filters wavelet , scale coefficients, add them together, , original list.

i can work haar wavelet filter, when try use cdf 9/7 filter doesn't produce same input. resulting list , original list sum same thing, however.

i'm sure it's stupid error in convolution, can't figure out. i've tried bunch of permutations of convolution, centering filter on index "i", instead of starting left edge of it, nothing has seemed work... it's 1 of bugs make me slap head when figure out.

here's code:

import random import math length = 128 array = list() row = list() scalecoefficients = list() waveletcoefficients = list() reconstruction = list()  def upsample(lst, index):     if (index % 2 == 0):         return 0.0     else:         return lst[index/2]  in range(length):     array.append(random.random())  ## cdf 9/7 wavelet (doesn't work?) dwtanalysislowpass = [.026749, -.016864, -.078223, .266864, .602949, .266864, -.078223, -.016864, .026749] in range(len(dwtanalysislowpass)):     dwtanalysislowpass[i] = math.sqrt(2.0) * dwtanalysislowpass[i] dwtanalysishighpass = [0.0, .091272, -.057544, -0.591272, 1.115087, -.591272, -.057544, .091272, 0.0] in range(len(dwtanalysishighpass)):     dwtanalysishighpass[i] = 1.0/math.sqrt(2.0) * dwtanalysishighpass[i]  dwtsynthesislowpass = [0.0, -.091272, -.057544, 0.591272, 1.115087, .591272, -.057544, -.091272, 0.0] in range(len(dwtsynthesislowpass)):     dwtsynthesislowpass[i] = 1.0/math.sqrt(2.0) * dwtsynthesislowpass[i] dwtsynthesishighpass = [.026749, .016864, -.078223, -.266864, .602949, -.266864, -.078223, .016864, .026749] in range(len(dwtsynthesishighpass)):     dwtsynthesishighpass[i] = math.sqrt(2.0) * dwtsynthesishighpass[i]  ## haar wavelet (works) ## c = 1.0/math.sqrt(2) ## dwtanalysislowpass = [c,c] ## dwtanalysishighpass = [c, -c] ## dwtsynthesislowpass = [c, c] ## dwtsynthesishighpass = [-c, c]   ## forward transform - need on half elements in range(0,length,2):     newval = 0.0     ## convolve next j elements     j in range(len(dwtanalysislowpass)):         index = + j         if(index >= length):             index = index - length          newval = newval + array[index]*dwtanalysislowpass[j]      scalecoefficients.append(newval)      newval = 0.0     j in range(len(dwtanalysishighpass)):         index = + j         if(index >= length):             index = index - length          newval = newval + array[index]*dwtanalysishighpass[j]      waveletcoefficients.append(newval)  ## inverse transform in range(length):     newval = 0.0     j in range(len(dwtsynthesishighpass)):         index = + j         if(index >= length):             index = index - length          newval = newval + upsample(waveletcoefficients, index)*dwtsynthesishighpass[j]      j in range(len(dwtsynthesislowpass)):         index = + j         if(index >= length):             index = index - length          newval = newval + upsample(scalecoefficients, index)*dwtsynthesislowpass[j]      reconstruction.append(newval)  print sum(reconstruction) print sum(array) print reconstruction print array 

incidentally, took filter values appendix here: http://www1.cs.columbia.edu/~rso2102/awr/files/overbeck2009awr.pdf, i've seen them used in bunch of matlab sample code well.

actually solved myself comparing coefficients, , reconstruction, code lifting implementation:

http://www.embl.de/~gpau/misc/dwt97.c

basically, 1) made boundary conditions symmetric, instead of periodic 2) had offset convolutions (and upsampling) in ways line up.

here's code in case else runs across problem. feel still over-complicating it, because it's not documented anywhere, @ least works. includes "switch" used test against reference, , had modify haar wavelet make work.

import random import math length = int() array = list() row = list() scalecoefficients = list() waveletcoefficients = list() reconstruction = list() switch = false  def upsample1(lst, index):     if (index % 2 == 0):         return lst[index/2]     else:         return 0.0  def upsample2(lst, index):     if (index % 2 == 0):         return 0.0     else:         return lst[index/2]  ## generate random list of floating point numbers if (not switch):     length = 128     in range(length):         array.append(random.random()) else:     length = 32     in range(32):         array.append(5.0+i+.4*i*i-.02*i*i*i)  ## first part calculates filters ## cdf 9/7 wavelet dwtanalysislowpass = [.026749, -.016864, -.078223, .266864, .602949, .266864, -.078223, -.016864, .026749] in range(len(dwtanalysislowpass)):     dwtanalysislowpass[i] = math.sqrt(2.0) * dwtanalysislowpass[i] dwtanalysishighpass = [.091272, -.057544, -0.591272, 1.115087, -.591272, -.057544, .091272] in range(len(dwtanalysishighpass)):     dwtanalysishighpass[i] = dwtanalysishighpass[i]/math.sqrt(2.0)  dwtsynthesislowpass = [-.091272, -.057544, 0.591272, 1.115087, .591272, -.057544, -.091272] in range(len(dwtsynthesislowpass)):     dwtsynthesislowpass[i] = dwtsynthesislowpass[i]/math.sqrt(2.0) dwtsynthesishighpass = [.026749, .016864, -.078223, -.266864, .602949, -.266864, -.078223, .016864, .026749] in range(len(dwtsynthesishighpass)):     dwtsynthesishighpass[i] = math.sqrt(2.0) * dwtsynthesishighpass[i]  ## haar wavelet ## c = 1.0/math.sqrt(2) ## dwtanalysislowpass = [c,c] ## dwtanalysishighpass = [c, -c] ## dwtsynthesislowpass = [-c, c] ## dwtsynthesishighpass = [c, c]  # forward transform. can skip every other sample since # removed in downsampling anyway in range(0,length,2):     newval = 0.0     ## convolve next j elements low-pass analysis filter     j in range(len(dwtanalysislowpass)):         index = + j - len(dwtanalysislowpass)/2         if(index >= length):             index = 2*length - index - 2         elif (index < 0):             index = -index          newval = newval + array[index]*dwtanalysislowpass[j]      # append new value list of scale coefficients     scalecoefficients.append(newval)      newval = 0.0     # convolve next j elements high-pass analysis filter     j in range(len(dwtanalysishighpass)):         index = + j - len(dwtanalysishighpass)/2 + 1         if(index >= length):             index = 2*length - index - 2         elif (index < 0):             index = -index          newval = newval + array[index]*dwtanalysishighpass[j]      # append new value list of wavelet coefficients     waveletcoefficients.append(newval)  # inverse transform in range(length):     newval = 0.0     # convolve upsampled wavelet coefficients high-pass synthesis filter     j in range(len(dwtsynthesishighpass)):         index = + j - len(dwtsynthesishighpass)/2         if(index >= length):             index = 2*length - index - 2         elif (index < 0):             index = -index          newval = newval + upsample2(waveletcoefficients, index)*dwtsynthesishighpass[j]      # convolve upsampled scale coefficients low-pass synthesis filter, ,     # add previous convolution     j in range(len(dwtsynthesislowpass)):         index = + j - len(dwtsynthesislowpass)/2         if(index >= length):             index = 2*length - index - 2         elif (index < 0):             index = -index          newval = newval + upsample1(scalecoefficients, index)*dwtsynthesislowpass[j]      reconstruction.append(newval)  print ("sums: ") print sum(reconstruction) print sum(array) print ("original signal: ") print array if (not switch):     print ("wavelet coefficients: ")     in range(len(scalecoefficients)):         print ("sc[" + str(i) + "]: " + str(scalecoefficients[i]))     in range(len(waveletcoefficients)):         print ("wc[" + str(i) + "]: " + str(waveletcoefficients[i])) print ("reconstruction: ") print reconstruction