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?
