CDF 9/7 Diskrete Wavelet-Transformation (Faltung)

Question

Ich versuche, ein einfaches eigenständiges Programm zu schreiben, das eine diskrete Wavelet-Transformationsebene auf einer 1D-Liste mit den CDF 9/7-Wavelets durchführt und dann rekonstruiert. Ich verwende nur die Faltungs/Filterbank-Methode, um zu verstehen, wie es funktioniert. Mit anderen Worten, falte die Liste mit einem Filter, um die Skalarkoeffizienten zu erhalten, falte die Liste mit einem anderen Filter, um die Wavelet-Koeffizienten zu erhalten, aber mache dies nur, indem du mit jedem zweiten Element beginnst. Dann vergrößere sie (d.h. füge Nullen zwischen die Elemente ein), wende Filter auf die Wavelet- und Skalarkoeffizienten an, addiere sie zusammen und erhalte die originale Liste.

Ich kann dies für den Haar-Wavelet-Filter zum Laufen bringen, aber wenn ich den CDF 9/7 Filter verwende, produziert er nicht die gleiche Eingabe. Die resultierende Liste und die ursprüngliche Liste summieren sich jedoch zu dem gleichen Wert.

Ich bin sicher, dass es ein sehr dummer Fehler bei der Faltung ist, aber ich kann einfach nicht darauf kommen. Ich habe viele Permutationen der Faltung ausprobiert, wie das Zentrieren des Filters auf dem Index "i", anstatt den linken Rand davon zu starten, aber nichts scheint zu funktionieren... Es ist wahrscheinlich einer dieser Fehler, die mich auf die Stirn schlagen werden, wenn ich ihn herausfinde.

Hier ist der 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 (funktioniert nicht?)
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 (funktioniert)
## c = 1.0/math.sqrt(2)
## DWTAnalysisLowpass = [c,c]
## DWTAnalysisHighpass = [c, -c]
## DWTSynthesisLowpass = [c, c]
## DWTSynthesisHighpass = [-c, c]

## Führe die Vorwärtstransformation durch - wir müssen dies nur für die Hälfte der Elemente tun
for i in range(0,length,2):
    newVal = 0.0
    ## Faltung der nächsten j-Elemente
    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)

## Führe die inverse Transformation durch
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

Übrigens habe ich die Filterwerte aus dem Anhang hier genommen: http://www1.cs.columbia.edu/~rso2102/AWR/Files/Overbeck2009AWR.pdf, aber ich habe sie auch in einer Menge Matlab-Beispielscode verwendet.

Answer 1

Actually I solved it myself by comparing the coefficients, and then the reconstruction, with the code from this lifting implementation:

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

Grundsätzlich habe ich 1) die Randbedingungen symmetrisch gemacht, anstatt periodisch 2) Ich musste die Faltung (und das Upsampling) auf bestimmte Weise verschieben, um alles auszurichten.

Hier ist der Code, für den Fall, dass jemand anderes auf das Problem stößt. Ich habe das Gefühl, dass dies immer noch zu kompliziert ist, insbesondere weil es nirgendwo wirklich dokumentiert ist, aber zumindest funktioniert es. Dies enthält auch den "Schalter", den ich zum Testen gegen diesen Verweis verwendet habe, und ich musste das Haar-Wavelet modifizieren, damit es funktioniert.

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]

## Generiere eine Zufallsliste von Gleitkommazahlen
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)

## Erster Teil berechnet nur die Filter
## 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]

# Führe die Vorwärtstransformation durch. Wir können jeden zweiten Wert überspringen, da er eh beim Downsampling entfernt wird
for i in range(0,length,2):
    newVal = 0.0
    ## Faltung der nächsten j Elemente durch den Tiefpass-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]

    # Hänge den neuen Wert an die Liste der Skalen-Koeffizienten an
    scaleCoefficients.append(newVal)

    newVal = 0.0
    # Faltung der nächsten j Elemente durch den Highpass-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]

    # Hänge den neuen Wert an die Liste der Wavelet-Koeffizienten an
    waveletCoefficients.append(newVal)

# Führe die Rücktransformation durch
for i in range(length):
    newVal = 0.0
    # Faltung der hinaufgesampelten Wavelet-Koeffizienten mit dem Highpass-Synthese-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]

    # Faltung der hinaufgesampelten Skalen-Koeffizienten mit dem Tiefpass-Synthese-Filter und
    # füge es der vorherigen Faltung hinzu
    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 ("Summen: ")
print sum(reconstruction)
print sum(array)
print ("Ursprüngliches Signal: ")
print array
if (not switch):
    print ("Wavelet-Koeffizienten: ")
    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 ("Rekonstruktion: ")
print reconstruction