How to implement a filter associated to a specific wavelet

Question

I am working on ECG signals, to eventually extract features in order to detect an arrhythmia and classify it. I am using Discrete Wavelet Transform with biorthogonal wavelet bior6.8 During my research, I came to know that wavelet transform is the convolution of the input processed signal with the daughter wavelets to get approximation and detail

$$X(a,b)=\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty} \Psi\left( \frac{t-b}{a}\right)x(t)dt $$ where $a$ is scaling and $b$ is time.

I can't find the expression of the mother wavelet anywhere also when it came to practice DWT is usually presented as a filter bank of high pass and low pass filter The question is what is the difference between the different wavelets if it is always presented as a bank of filters

In my work, I used these two Butterworth high pass and low pass filters, but I still can't explain my choice, I read that Butterworth is the most used in signal processing and that it optimizes the frequency response in the passband, getting as much as you can from the wanted frequency Still, I have no arguments why I shouldn't use any others and not sure whether it is the correct way to implement bior 6.8 wavelet as I do not know any other way to implement wavelets and I would implement Daubechies or any other the same which do not make sense

from scipy.signal import filtfilt , butter

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data,padlen=0)
    return y

---

from scipy.signal import butter, lfilter

def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

So basically and to resume, the question is what is the difference between wavelet transforms and how can we value that during implementation and then, the use of BF in this type of wavelet is it correct?

NB: the image inserted is from Wikipedia

Answer 1

The integral formula is the continuous wavelet transform, where you indeed convolve the data with a (mother) wavelet at all locations and scales. It is, by essence, redundant, and very mild conditions are required from continuous wavelet shapes to have both analysis and synthesis, and many shapes exist in analytic form.

If you want to discretize the above scheme in a critical way, ie to convert $N$ samples of a signal into $N$ coefficients, then you reach the realm of the discrete wavelet transform. Here, quite a limited choice of admissible wavelets exists (compared to the continuous case). And there, most continuous wavelets don't fit in the frame (pun intended). In other words, they cannot be discretized perfectly. The implementation technique was reversed: it was shown that the DWT could be implemented as the filter bank you mention. And from wavelet properties, filters properties where designed. Then, given a specific pair of admissible filters, wavelets are obtained through an iterative process, and generally, there is not closed form or formulae for such wavelets. If you are interested in more details, you can read Filters and Wavelets for Dyadic Analysis by D. S. G. Pollock:

Except in the polar cases of the Haar wavelet and the Shannon or sinc function wavelet, an explicit functional form for the wavelet or the scaling function is unlikely to be available. Nor is there, in most practical applications, an indispensable requirement to represent of these functions graphically

A whole DWT framework exists in Python: PyWavelets - Wavelet Transforms in Python. For instance, you can get filter coefficients for the biorthogonal 6.8 wavelet [EDIT: on the page wavelets.pybytes.com/wavelet/bior6.8 under Coefficients, click "Show values" (or "Copy")].

Conversely, if you take any pair of low/high pass filters, it will not necessarily become a genuine DWT filter bank, because they don't have the necessary "wavelet type" properties.

So, basically, to implement a DWT the standard way, you should pick filter coefficients tabulated in the literature, or listed at the PyWavelet property browser.

Answer 2

You've mentioned Butterworth filters for doing the wavelet analysis using bior6.8. If you want to perform the Discrete Wavelet Transform using some specific wavelet, then you must use its Perfect Reconstruction Filter Bank. Each wavelet function has its associated set of filter values for decomposition and reconstruction - they are calculated from the Mother and Daughter wavelet.

Since you were using Python, I am going to give you an example using PyWavelets. It is as simple as creating a wavelet object and asking for the filter bank. Here is a plot of the filter values:

Since these are FIR filter coefficients, then we can calculate the amplitude response by taking the DFT of the coefficient values:

---

Values of filter coefficients:

dec_lo = [0.0,
 0.0019088317364812906,
 -0.0019142861290887667,
 -0.016990639867602342,
 0.01193456527972926,
 0.04973290349094079,
 -0.07726317316720414,
 -0.09405920349573646,
 0.4207962846098268,
 0.8259229974584023,
 0.4207962846098268,
 -0.09405920349573646,
 -0.07726317316720414,
 0.04973290349094079,
 0.01193456527972926,
 -0.016990639867602342,
 -0.0019142861290887667,
 0.0019088317364812906]

dec_hi = [-0.0,
 0.0,
 -0.0,
 0.014426282505624435,
 -0.014467504896790148,
 -0.07872200106262882,
 0.04036797903033992,
 0.41784910915027457,
 -0.7589077294536541,
 0.41784910915027457,
 0.04036797903033992,
 -0.07872200106262882,
 -0.014467504896790148,
 0.014426282505624435,
 -0.0,
 0.0,
 -0.0,
 0.0]

rec_lo = [0.0,
 0.0,
 0.0,
 0.014426282505624435,
 0.014467504896790148,
 -0.07872200106262882,
 -0.04036797903033992,
 0.41784910915027457,
 0.7589077294536541,
 0.41784910915027457,
 -0.04036797903033992,
 -0.07872200106262882,
 0.014467504896790148,
 0.014426282505624435,
 0.0,
 0.0,
 0.0,
 0.0]

rec_hi = [0.0,
 -0.0019088317364812906,
 -0.0019142861290887667,
 0.016990639867602342,
 0.01193456527972926,
 -0.04973290349094079,
 -0.07726317316720414,
 0.09405920349573646,
 0.4207962846098268,
 -0.8259229974584023,
 0.4207962846098268,
 0.09405920349573646,
 -0.07726317316720414,
 -0.04973290349094079,
 0.01193456527972926,
 0.016990639867602342,
 -0.0019142861290887667,
 -0.0019088317364812906]

---

A complementary code which might be useful for you to investigate other wavelets:

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

# Create wavelet and extract the filters
wavelet_name = 'bior6.8'
wavelet = pywt.Wavelet(wavelet_name)
dec_lo, dec_hi, rec_lo, rec_hi = wavelet.filter_bank

# Filter coefficients
plt.figure()
plt.subplot(221)
plt.stem(dec_lo)
plt.grid()
plt.title('{} low-pass decomposition filter'.format(wavelet_name))
plt.subplot(222)
plt.stem(dec_hi)
plt.grid()
plt.title('{} high-pass decomposition filter'.format(wavelet_name))
plt.subplot(223)
plt.stem(rec_lo)
plt.grid()
plt.title('{} low-pass reconstruction filter'.format(wavelet_name))
plt.subplot(224)
plt.stem(rec_hi)
plt.grid()
plt.title('{} high-pass reconstruction filter'.format(wavelet_name))

# Frequency responses
dec_lo_fr = np.abs(np.fft.rfft(dec_lo, 128))
dec_hi_fr = np.abs(np.fft.rfft(dec_hi, 128))
rec_lo_fr = np.abs(np.fft.rfft(rec_lo, 128))
rec_hi_fr = np.abs(np.fft.rfft(rec_hi, 128))

plt.figure()
plt.subplot(211)
plt.plot(dec_lo_fr, label='Low-pass')
plt.hold(True)
plt.plot(dec_hi_fr, label='High-pass')
plt.grid()
plt.legend()
plt.title('Frequency responses of {} decomposition filters'.format(wavelet_name))
plt.subplot(212)
plt.plot(rec_lo_fr, label='Low-pass')
plt.hold(True)
plt.plot(rec_hi_fr, label='High-pass')
plt.grid()
plt.legend()
plt.title('Frequency responses of {} reconstruction filters'.format(wavelet_name))

plt.show()

---

Regarding the actual filtering process you can either use scipy.signal.lfilter or np.convolve (x being the input signal):

y = lfilter(dec_lo, 1, x)

or

y = convolve(x, dec_lo)

I would suggest using PyWavelets to do all of that for you, instead of re-inventing the wheel. However, If you really want to implement a single level of DWT on your own then I second the numpy option. Here is a complete example:

t = np.linspace(0, 1.0, 128)
x = np.sin(2*np.pi*10*t)

# Perform manual DWT and decimate
cA = np.convolve(x, dec_lo)[1::2]
cD = np.convolve(x, dec_hi)[1::2]

plt.figure()
plt.subplot(211)
plt.plot(cA)
plt.grid()
plt.title('Approximation coefficients')

plt.subplot(212)
plt.plot(cD)
plt.grid()
plt.title('Detail coefficients')

plt.show()

Which gives: