Feature extraction/reduction using DWT

Question

For a given time series which is n timestamps in length, we can take Discrete Wavelet Transform (using 'Haar' wavelets), then we get (for an example, in Python) -

>>> import pywt
>>> ts = [2, 56, 3, 22, 3, 4, 56, 7, 8, 9, 44, 23, 1, 4, 6, 2]
>>> (ca, cd) = pywt.dwt(ts,'haar')
>>> ca
array([ 41.01219331,  17.67766953,   4.94974747,  44.54772721,
        12.02081528,  47.37615434,   3.53553391,   5.65685425])
>>> cd
array([-38.18376618, -13.43502884,  -0.70710678,  34.64823228,
        -0.70710678,  14.8492424 ,  -2.12132034,   2.82842712])

where ca and cd are approximation and detailed coefficients. Now if I use all of them I can construct my original time series back using inverse DWT. But instead I want to use a fewer coefficients (like in Fourier Transform if we use only first few coefficients, we can approximately reconstruct the original time series). If I just use ca or just use cd I don't get the desired results. If I use only we coefficients from each of them (like first 4), I get only half of the time series.

How should I select the coefficients (from ca and cd) so that I can approximately create the original signal from them (i.e. most of its energy)?

Answer 1

I think it is kind'a similar to soft and hard thresholding using in wavelet de-noising. Have you come across this topic? pywt has already an in-built function for this purpose. Please take a closer look at this code and try to play with it:

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

ts = [2, 56, 3, 22, 3, 4, 56, 7, 8, 9, 44, 23, 1, 4, 6, 2]

(ca, cd) = pywt.dwt(ts,'haar')

cat = pywt.thresholding.soft(ca, np.std(ca)/2)
cdt = pywt.thresholding.soft(cd, np.std(cd)/2)

ts_rec = pywt.idwt(cat, cdt, 'haar')

plt.close('all')

plt.subplot(211)
# Original coefficients
plt.plot(ca, '--*b')
plt.plot(cd, '--*r')
# Thresholded coefficients
plt.plot(cat, '--*c')
plt.plot(cdt, '--*m')
plt.legend(['ca','cd','ca_thresh', 'cd_thresh'], loc=0)
plt.grid('on')

plt.subplot(212)
plt.plot(ts)
plt.hold('on')
plt.plot(ts_rec, 'r')
plt.legend(['original signal', 'reconstructed signal'])
plt.grid('on')
plt.show()

This will produce following - hope this will guide you.

Note on running the code with the latest versions:

if you get an error like this:

module 'pywt' has no attribute 'thresholding'

use on line 9 and 10:

 cat = pywt.threshold(ca, np.std(ca)/2, mode='soft')
 cdt = pywt.threshold(cd, np.std(cd)/2, mode='soft')

if you get an error on:

plt.hold('on') then comment this line out:

  # plt.hold('on')

Answer 2

A notion that fits wavelets well is that of NLA, non-linear approximation. Given a length $N$ signal $[x_n]$, and its transformation coefficients $[X_m]$ (of length $M$) under transform $\mathbb{T}$. The best $K$-term approximation would be a subset of $K$ terms of indices from $[X_m]$, denoted by $\sigma(k)$: $[X_\sigma(k)]$, $1\le k\le K$. It will be best under some metric $\mu$, measuring the difference between $x$ and $\hat{x}^\mathbb{T}_K$, the recovered signal from the "candidate" $K$-term approximation, all non-kept coefficients being set to zero:

$$\mu (x-\hat{x}^\mathbb{T}_K)$$

What is important, and was the subject of many reasearch, is how good are those approximations (how fast its decays with $K$)

• given a class of signals (e.g. piecewise polynomial)
• given a family of transformations.

The overall problem is very complicated, but let us made a long story short. A nice heuristic for local energy-preserving transformations (like wavelets) is to pick the highest coefficients in magnitude. In other words, the $[X_\sigma(k)]$, $1\le k\le K$ are the largest coefficients. This is consistent with the answer by jojek; all coefficients above a threshold are somehow the largest. In 1D, this is neat because wavelets both "compact" smooth/regular parts in the cA coefficients, and non-smoothness is well-concentrated (sparse) in cD coefficients.

In practice, it is well documented that it can be important to normalize coefficients across scales. Or to select coefficients jointly: not independently in scale and location, only looking at magnitude. This can be called non-scalar thresholding, group thresholding, etc.

All of the above is good all long as large coefficients well approximate, or suit your application. When doing feature extraction, it might be useful to first identify, or learn, what coefficients/bands of your wavelet transform are indeed useful to you. Two proposed steps:

1. with proper coefficient normalization (if needed), verify if picking the highest coefficients is efficient for your purpose
2. if not, use the physics to select educated guesses, or learn features according to your target.