I'm trying to write a simple self-contained program that does a single level of a discrete wavelet transform on a 1D list, using the CDF 9/7 wavelets, and then reconstructs it. I'm just using the convolution/filter-bank method to get a grasp on how it works. In other words, convolve the list with a filter to get the scale coefficients, convolve the list with a different filter to get the wavelet coefficients, but only do this starting with every other element. Then upsample (i.e. add zeroes in between the elements), apply filters to the wavelet and scale coefficients, add them together, and get the original list.
I can get this to work for the Haar wavelet filter, but when I try to use the CDF 9/7 filter it doesn't produce the same input. The resulting list and the original list do sum to the same thing, however.
I'm sure it's a very stupid error in the convolution, but I just can't figure it out. I've tried a bunch of permutations of the convolution, like centering the filter on the index "i", instead of starting the left edge of it, but nothing has seemed to work... It's probably one of those bugs that will make me slap my head when I figure it out.
Here's the 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]
for i in range(length):
array.append(random.random())
## CDF 9/7 Wavelet (doesn't work?)
DWTAnalysisLowpass = [.026749, -.016864, -.078223, .266864, .602949, .266864, -.078223, -.016864, .026749]
for i 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]
for i 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]
for i in range(len(DWTSynthesisLowpass)):
DWTSynthesisLowpass[i] = 1.0/math.sqrt(2.0) * DWTSynthesisLowpass[i]
DWTSynthesisHighpass = [.026749, .016864, -.078223, -.266864, .602949, -.266864, -.078223, .016864, .026749]
for i 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]
## Do the forward transform - we only need to do it on half the elements
for i in range(0,length,2):
newVal = 0.0
## Convolve the next j elements
for j in range(len(DWTAnalysisLowpass)):
index = i + j
if(index >= length):
index = index - length
newVal = newVal + array[index]*DWTAnalysisLowpass[j]
scaleCoefficients.append(newVal)
newVal = 0.0
for j in range(len(DWTAnalysisHighpass)):
index = i + j
if(index >= length):
index = index - length
newVal = newVal + array[index]*DWTAnalysisHighpass[j]
waveletCoefficients.append(newVal)
## Do the inverse transform
for i in range(length):
newVal = 0.0
for j in range(len(DWTSynthesisHighpass)):
index = i + j
if(index >= length):
index = index - length
newVal = newVal + upsample(waveletCoefficients, index)*DWTSynthesisHighpass[j]
for j in range(len(DWTSynthesisLowpass)):
index = i + 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, I took the filter values from the appendix here: http://www1.cs.columbia.edu/~rso2102/AWR/Files/Overbeck2009AWR.pdf, but I've seen them used in a bunch of matlab sample code as well.
numpy
so you can do things likeDWTSynthesisHighpass *= math.sqrt(2.0)
w/o a loop. More generally, you almost never need to dofor i in range(len(stuff))
when you can dofor x in stuff
. – Andrew Jaffefor x in stuff
; I could use slicing in the actual transforms but I think that would obfuscate it more. – Andrew