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)) x_double = np.linspace(x, x[-1], x.shape*2) f = interpolate.interp1d(x, data, kind='linear') data_upsampled = f(x_double) return data_upsampled # Upsampling by interpolation, reconstruction cA = coefs_own 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.plot(coefs_own) axes.plot(coefs_own) axes.set_title('Own') axes.legend(['Approx.', 'Detail']) axes.plot(coefs) axes.plot(coefs) axes.set_title('PyWavelets') axes.legend(['Approx.', 'Detail']) # Plot both reconstructions fig, axes = plt.subplots(3, 1) axes.plot(sig) axes.set_title('Original') axes.plot(sig_rec_own) axes.set_title('Self-made reconstruction') axes.plot(sig_rec_pywt) axes.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?