Discrete Wavelet Transform - Visualizing Relation between Decomposed Detail Coefficients and Signal

Question

I'm trying to directly visualize the relation between discrete wavelet transform (DWT) detail coefficients and the original signal/its reconstruction. The goal is to show their relation in an intuitive way. I would like to ask (see questions below): if the idea and process I've come up with is correct so far, and if I am right that it might be better subtract the 1st level approximation from the original signal before visualizing their relation.

Minimal example

Here's the minimal example I base my explanation on, using the ECG example data of Python's pywavelets, which has 1024 values, as a simple 1D signal:

import pywt
import pywt.data
import numpy as np
import matplotlib.pyplot as plt

x = pywt.data.ecg()
plt.plot(x)
plt.legend(['Original signal'])

Decomposition is done using a Symmlet 5 with a total of 6 levels:

w = pywt.Wavelet('sym5')
plt.plot(w.dec_lo)
coeffs = pywt.wavedec(x, w, level=6)

(Lossy) reconstruction of the signal works as expected when intentionally leaving out detail coefficients of higher levels (signals are plotted on uniform x-scale [0,1] for convenience):

def reconstruction_plot(yyy, **kwargs):
    """Plot signal vector on x [0,1] independently of amount of values it contains."""
    plt.plot(np.linspace(0, 1, len(yyy)), yyy, **kwargs)

reconstruction_plot(pywt.waverec(coeffs, w)) # full reconstruction 
#reconstruction_plot(pywt.waverec(coeffs[:-1] + [None] * 1, w)) # leaving out detail coefficients up to lvl 5
#reconstruction_plot(pywt.waverec(coeffs[:-2] + [None] * 2, w)) # leaving out detail coefficients up to lvl 4
#reconstruction_plot(pywt.waverec(coeffs[:-3] + [None] * 3, w)) # leaving out detail coefficients up to lvl 3
reconstruction_plot(pywt.waverec(coeffs[:-4] + [None] * 4, w)) # leaving out detail coefficients up to lvl 2
#reconstruction_plot(pywt.waverec(coeffs[:-5] + [None] * 5, w)) # leaving out detail coefficients up to lvl 1
reconstruction_plot(pywt.waverec(coeffs[:-6] + [None] * 6, w)) # leaving out all detail coefficients = reconstruction using lvl1 approximation only
plt.legend(['Full reconstruction', 'Reconstruction using detail coefficients lvl 1+2', 'Reconstruction using lvl 1 approximation only'])

The DWT above yields a level 1 approximation vector of 24 values, level 1 details coefficient vector of 24 values, level 2 detail vector of 40 values, level 3 of 72 values, level 4 of 135 values, level 5 of 262 values, and level 6 of 516 values:

plt.stem(coeffs[1]); plt.legend(['Lvl 1 detail coefficients'])
plt.stem(coeffs[2]); plt.legend(['Lvl 2 detail coefficients'])
plt.stem(coeffs[3]); plt.legend(['Lvl 3 detail coefficients'])
plt.stem(coeffs[4]); plt.legend(['Lvl 4 detail coefficients'])
plt.stem(coeffs[5]); plt.legend(['Lvl 5 detail coefficients'])
plt.stem(coeffs[6]); plt.legend(['Lvl 6 detail coefficients'])

It seems we see clear patterns around the spikes in original signals (also pay attention to the y scale of the plots above).

Now to my questions:

1. It is correct that we can directly relate those coefficients to the signal? Amplitude of the coefficient corresponds to amplitude with which the wavelet occurs in the signal (y axis), and position of the coefficient corresponds to the time (x axis). Or is there something in between we need to consider?

2. After the DWT the final lvl1 approximation remains. Does it make sense to not visualize the relation of the detail coefficients with the original signal, but instead with the original signal minus the lvl1 approximation? (I know that I would most likely also see the relation between coefficients and signal without doing this, see e.g. plots below. It's just for it this makes sense or not. Should it makes sense for lvl1 detail coefficients then it might also make sense for lvl2 detail coefficients to be compared to the original signal minus the lvl2 approximation, right?). An example:

   # Reconstruction of signal using just lvl1 approximation
   approx_lvl1 = pywt.waverec(coeffs[:-6] + [None] * 6, w)
   # interpolate to original amount of samples (necessary due to numeric solution of transformation not yielding same amount of values)
   approx_lvl1_interp = np.interp(x=np.arange(0, 1024), xp=np.linspace(0, 1024, len(approx_lvl1)), fp=approx_lvl1)
   x_without_lvl1approx = x - approx_lvl1_interp
   

3. The direct visualization of the relation between detail coefficients and signal I use just plots both the signal and the coefficients on an x axis of [0,1]. This should conceptually be valid, but I am unsure if I would actually need an offset towards the margins (e.g. first and last coefficient of the vector not being positioned at the very beginning or end of the signal):

   def reconstruction_stem(yyy, **kwargs):
       """Plot coefficient vector on x [0,1] independently of amount of values it contains."""
       plt.stem(np.linspace(0, 1, len(yyy)), yyy, **kwargs)
   
   reconstruction_plot(x, color='orange')
   reconstruction_plot(x_without_lvl1approx, color='red')
   reconstruction_stem(coeffs[1])
   plt.legend(['Original signal', 'Original signal - lvl1 approximation', 'Detail coefficients'])
   

Is there an intuitive explanation for the strong coefficients that are not directly at the positions of peaks in the original data (e.g. at level 1 the lowest one (strongest negative) at around 0.25, as well as the highest one (strongest positive) around 0.75)? Though there is a clear pattern (positive lag+negative amplitude, negative lag+positive amplitude) those seem a bit "far off" for me. But there is probably a good explanation for that.

Thanks for answering!

Answer 1

It is necessary to distinguish clearly between the approximation and detail coefficients at each decomposition level and the levels associated Details and Approximation which involves not only the coefficients but also the inverse filters at that level respectively

Answer 2

I am only now starting to dabble with wavelets, and am still struggling even with very basic questions like "how does one choose from the gaggle of available wavelets" (probably has to do with the number of levels you need to achieve "good enough" representation), and "what is all the hoopla about denoising with wavelets", because I seem to be able to achieve better results for my type of data with straight gaussian denoising or median filters. But I digress....

One thing I noted in the above is that your level numbering seems inconsistent with what I believe I the usual wavelet conventions. In particular coeffs[0] is the approximation amplitudes at the last level, in your case 6 coeffs[1] is the detail amplitude at level 6 coeffs[2] is the detail amplitude at level 5 ... coeffs[6] is the detail amplitude at level 1

So your reconstructions are from level 5 and 6 only, not level 1 and 2 as indicated in your plots.

=========

Update: I messed with your code more, and I think your idea of illustrating the correlation between coefficients and signal features is sound, but not perfect. I have mucked around with your code a bit to illustrate this better, see below. Note that at each step I rescale the coefficients to the magnitude of the signal. This allows one to talk about the concept of thresholding too.

import pywt
import pywt.data
import numpy as np
import matplotlib.pyplot as plt

plt.close('all')

def reconstruction_plot(yyy, **kwargs):
    """Plot signal vector on x [0,1] independently of amount of values it contains."""
    #plt.figure()
    #plt.plot(np.linspace(0, 1, len(yyy)), yyy, **kwargs)
    ym = np.median(yyy)
    plt.plot(np.linspace(0, 1., num=len(yyy)), yyy-ym, **kwargs)


def reconstruction_stem(yyy, xmax, **kwargs):
    """Plot coefficient vector on x [0,1] independently of amount of values it contains."""
    ymax = yyy.max()
    plt.stem(np.linspace(0, 1., num=len(yyy)), yyy*(xmax/ymax), **kwargs)


x = pywt.data.ecg()
w = pywt.Wavelet('sym5')
nl = 6
coeffs = pywt.wavedec(x, w, level=nl)


'''
plt.figure()
plt.stem(coeffs[1]); plt.legend(['Lvl 6 detail coefficients'])
plt.figure()
plt.stem(coeffs[2]); plt.legend(['Lvl 5 detail coefficients'])
plt.figure()
plt.stem(coeffs[3]); plt.legend(['Lvl 4 detail coefficients'])
plt.figure()
plt.stem(coeffs[4]); plt.legend(['Lvl 3 detail coefficients'])
plt.figure()
plt.stem(coeffs[5]); plt.legend(['Lvl 2 detail coefficients'])
plt.figure()
plt.stem(coeffs[6]); plt.legend(['Lvl 1 detail coefficients'])
'''


xmax = x.max()
for i in range(nl):
    plt.figure()
    reconstruction_plot(x) # original signal 
    #reconstruction_plot(pywt.waverec(coeffs, w)) # full reconstruction 
    reconstruction_plot(pywt.waverec(coeffs[:i+2] + [None] * (nl-i-1), w)) # partial reconstruction 
    reconstruction_stem(coeffs[i+1], xmax, markerfmt ='none', linefmt='r-')
    #plt.legend(['Original', 'Full reconstruction', ('Rec to lvl %d')%(nl-i), ('Details for lvl %d')%(nl-i)])
    plt.legend(['Original', ('Rec to lvl %d')%(nl-i), ('Details for lvl %d')%(nl-i)])