我试图编写一个简单的自包含程序,使用cdf9/7小波对一维列表进行单级离散小波变换,然后重建它。我只是用卷积/滤波器组方法来了解它的工作原理。换言之,用滤波器卷积列表以获得尺度系数,用不同的滤波器卷积列表以获得小波系数,但仅从其他元素开始。然后向上采样(即在元素之间加零),对小波系数和尺度系数应用滤波器,将它们相加,得到原始列表。在
我可以让Haar小波滤波器工作,但是当我尝试使用cdf9/7滤波器时,它不会产生相同的输入。但是,结果列表和原始列表的总和是相同的。在
我敢肯定这是卷积中一个非常愚蠢的错误,但我就是搞不懂。我尝试过卷积的一系列排列,比如把过滤器放在索引“I”上,而不是从它的左边缘开始,但似乎没有任何效果。。。这可能是其中一个虫子,当我弄清楚的时候,会让我拍拍我的头。在
代码如下:
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
顺便说一下,我从这里的附录中获取了过滤器值:http://www1.cs.columbia.edu/~rso2102/AWR/Files/Overbeck2009AWR.pdf,但我也看到它们在一堆matlab示例代码中使用过。在
实际上,我自己解决了这个问题,通过比较系数,然后重建,与这个提升实现的代码:
http://www.embl.de/~gpau/misc/dwt97.c
基本上,我 1) 使边界条件对称,而不是周期性的 2) 必须以一定的方式抵消卷积(和上采样)以使其对齐。在
下面是代码,以防其他人遇到问题。我觉得这仍然过于复杂化了,特别是因为它在任何地方都没有真正的文档记录,但至少它是有效的。这也包括我用来对照那个参考进行测试的“开关”,我必须修改Haar小波使其工作。在
相关问题 更多 >
编程相关推荐