0
$\begingroup$

Trying to decompose and reconstruct a signal using a to some extent self-made implementation of DWT for some reason fails. The result looks highpass filtered and/or shifted. I wanted to write the code in Python not using the PyWavelets package (except for getting the wavelet object).

import numpy as np
from scipy import interpolate, signal
import pywt
import matplotlib.pyplot as plt

# Test signal and wavelet
sig = np.zeros((10*512))
for i in range(256):
    sig += i * np.cos(2*np.pi*i*np.linspace(0, 10, 10*512)+i**2*2*np.pi/((i+1)/1024))
w = pywt.DiscreteContinuousWavelet('db6')

I convolved the test signal with the decomposition filters from the wavelet object and used scipy's signal.decimate(...). As far as I understood, the decomposition filters do not lowpass filter the signal (or do they and is this where I went wrong?). However, a simple downsampling instead of decimation did not yield better results.

# Downsampling by decimation, decomposition
coefs_own = []
cA = np.copy(sig)
for i in range(1):
    cA_last = cA
    cA = np.convolve(cA_last, w.dec_lo, 'same')
    cD = np.convolve(cA_last, w.dec_hi, 'same')
    cA = signal.decimate(cA, 2)
    cD = signal.decimate(cD, 2)
    coefs_own.insert(0, cD)
coefs_own.insert(0, cA)

And this is my reconstruction using linear interpolation and a convolution with the recomposition filters provided by the wavelet object.

# Helper function
def upsample(data):
    x = np.array(range(data.shape[0]))
    x_double = np.linspace(x[0], x[-1], x.shape[0]*2)
    f = interpolate.interp1d(x, data, kind='linear')
    data_upsampled = f(x_double)
    return data_upsampled

# Upsampling by interpolation, reconstruction
cA = coefs_own[0]
for i in range(1):
    cA_up = upsample(cA)
    cD_up = upsample(coefs_own[i+1])
    cA = np.convolve(cA_up, w.rec_lo, 'same') + np.convolve(cD_up, w.rec_hi, 'same')
sig_rec_own = cA

For comparison, see the multilevel decomposition by PyWavelets, which does a perfect reconstruction.

# PyWavelets
coefs = pywt.wavedec(sig, w, mode='zero', level=1)
sig_rec_pywt = pywt.waverec(coefs, w, mode='zero')

# Plot both decomposition versions
fig, axes = plt.subplots(2, 1)
axes[0].plot(coefs_own[0])
axes[0].plot(coefs_own[1])
axes[0].set_title('Own')
axes[0].legend(['Approx.', 'Detail'])
axes[1].plot(coefs[0])
axes[1].plot(coefs[1])
axes[1].set_title('PyWavelets')
axes[1].legend(['Approx.', 'Detail'])

# Plot both reconstructions
fig, axes = plt.subplots(3, 1)
axes[0].plot(sig)
axes[0].set_title('Original')
axes[1].plot(sig_rec_own)
axes[1].set_title('Self-made reconstruction')
axes[2].plot(sig_rec_pywt)
axes[2].set_title('PyWavelets reconstruction')

plt.show()

Does anybody have an idea, why the decomposition and the reconstruction both obviously are erroneous? What did I miss in my implementation?

enter image description here enter image description here

$\endgroup$
  • $\begingroup$ Why the linear interpolation? Usually the wavelet filters are designed that you fill with zeros while upsampling. $\endgroup$ – LutzL Apr 29 at 15:59
  • $\begingroup$ Thanks, you are right. It must be interpolated by inserting zeros. I would mark your answer as accepted, but I cannot find any button for that. As you see, I am newly registered here so I might have overlooked the accept-button? $\endgroup$ – Sergeant Salty May 2 at 13:42
  • $\begingroup$ I put a longer version of the comment as answer, as of now you can only vote on answers, later you can give kudos to comments. $\endgroup$ – LutzL May 2 at 13:55
0
$\begingroup$

The wavelet transform is usually constructed in such a way that in the reconstruction the input signals are either convoluted with each of the downsampled filters and then summed up, or equivalently, are upsampled filling the gaps with zeros.

Upsampling with interpolation is equal to applying an additional filter $[0.5,1,0.5]$ which is a low-pass, smoothening filter. Esp. on the high-frequency or detail component this will lead to large errors.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.