logo标签列表

CDF 9/7离散小波变换（卷积）

pythonwavelethaar-wavelet
3
我正在尝试编写一个简单的、自包含的程序，对一个一维列表进行一级离散小波变换，使用CDF 9/7小波，然后进行重构。我只是使用卷积/滤波器组方法来理解它的工作原理。换句话说，用一个滤波器对列表进行卷积以获得比例系数，用另一个滤波器对列表进行卷积以获得小波系数，但仅从每个其他元素开始执行此操作。然后上采样（即在元素之间添加零），将小波和比例系数应用于滤波器，将它们加在一起，得到原始列表。
我可以使用Haar小波滤波器使其正常工作，但当我尝试使用CDF 9/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 示例代码中使用。

2
提示：使用 numpy，这样您就可以像 DWTSynthesisHighpass *= math.sqrt(2.0) 这样做而不需要循环。更一般地说，当您可以使用 for x in stuff 时，您几乎永远不需要执行 for i in range(len(stuff))。 - Andrew Jaffe
是的，我知道，但是我当前的机器上没有安装numpy，所以想快速草拟一下。在这种情况下，我不能真正做到for x in stuff；我可以在实际转换中使用切片，但我认为那会更加混淆。 - Andrew
转念一想，如果我没有NumPy，使用map和lambda函数可能会更加简洁，但我认为生成过滤器并没有出现问题。 - Andrew
1个回答
2
实际上，我通过比较系数和重构与这个lifting实现的代码来解决了它。

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

基本上，我： 1）将边界条件改为对称，而不是周期性的 2）必须以某种方式偏移卷积（和上采样），以使其全部对齐。

以下是代码，以防其他人遇到此问题。我觉得这仍然过于复杂，尤其是因为它在任何地方都没有真正记录，但至少它有效。这也包括我用来测试参考文献的“开关”，我还必须修改Haar小波才能使其工作。

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 a random list of floating point numbers
if (not switch):
    length = 128
    for i in range(length):
        array.append(random.random())
else:
    length = 32
    for i in range(32):
        array.append(5.0+i+.4*i*i-.02*i*i*i)

## First Part Just Calculates the Filters
## CDF 9/7 Wavelet
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 = [.091272, -.057544, -0.591272, 1.115087, -.591272, -.057544, .091272]
for i in range(len(DWTAnalysisHighpass)):
    DWTAnalysisHighpass[i] = DWTAnalysisHighpass[i]/math.sqrt(2.0)

DWTSynthesisLowpass = [-.091272, -.057544, 0.591272, 1.115087, .591272, -.057544, -.091272]
for i in range(len(DWTSynthesisLowpass)):
    DWTSynthesisLowpass[i] = DWTSynthesisLowpass[i]/math.sqrt(2.0)
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
## c = 1.0/math.sqrt(2)
## DWTAnalysisLowpass = [c,c]
## DWTAnalysisHighpass = [c, -c]
## DWTSynthesisLowpass = [-c, c]
## DWTSynthesisHighpass = [c, c]

# Do the forward transform. We can skip every other sample since they would
# be removed in the downsampling anyway
for i in range(0,length,2):
    newVal = 0.0
    ## Convolve the next j elements by the low-pass analysis filter
    for j in range(len(DWTAnalysisLowpass)):
        index = i + j - len(DWTAnalysisLowpass)/2
        if(index >= length):
            index = 2*length - index - 2
        elif (index < 0):
            index = -index

        newVal = newVal + array[index]*DWTAnalysisLowpass[j]

    # append the new value to the list of scale coefficients
    scaleCoefficients.append(newVal)

    newVal = 0.0
    # Convolve the next j elements by the high-pass analysis filter
    for j in range(len(DWTAnalysisHighpass)):
        index = i + 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 the new value to the list of wavelet coefficients
    waveletCoefficients.append(newVal)

# Do the inverse transform
for i in range(length):
    newVal = 0.0
    # convolve the upsampled wavelet coefficients with the high-pass synthesis filter
    for j in range(len(DWTSynthesisHighpass)):
        index = i + 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 the upsampled scale coefficients with the low-pass synthesis filter, and
    # add it to the previous convolution
    for j in range(len(DWTSynthesisLowpass)):
        index = i + 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: ")
    for i in range(len(scaleCoefficients)):
        print ("sc[" + str(i) + "]: " + str(scaleCoefficients[i]))
    for i in range(len(waveletCoefficients)):
        print ("wc[" + str(i) + "]: " + str(waveletCoefficients[i]))
print ("Reconstruction: ")
print reconstruction
我使用不同的过滤器（例如Haar）运行了这段代码，但它对它们不起作用。它是否适用于任何类型的过滤器？ - sinner
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接