pythonsignal-processingfrequencywavelet-transform

How to combine Wavelet Transform and Frequency Filtering


I need to implement the following de-noising on ECG signal:

I do not know how to perform the second step in Python (PyWavelets), because I can modify only the detail and approximation coefficients and I do not know how to relate them to the frequencies.

How should I proceed?

This is my code

    import pywt

    #DWT
    coeff = pywt.wavedec(data,'db6',level=9)

    #filter the 0-0.35Hz frequencies in the 9-th level?


    #reconstruct the signal
    y = pywt.waverec( coeff[:8]+ [None] * 2, 'db6' )

Solution

  • My previous (now deleted) answer was a little confusing. Here I will try providing you with a hands-on example showing that reconstructing ECG data sampled at 360Hz using only the 'db6' approximation coefficients is (roughly) equivalent to low-pass filtering these data using a cut-off frequency of 0.35Hz.

    In the code example below, I import an ECG time series from scipy (from scipy.misc import electrocardiogram); they are sampled at 360Hz, just like yours. I will filter these data data using:

    Here's the code example:

    import pywt
    import numpy as np
    from scipy.misc import electrocardiogram
    import scipy.signal as signal
    import matplotlib.pyplot as plt
    
    wavelet_type='db6'
    data = electrocardiogram()
    
    DWTcoeffs = pywt.wavedec(data,wavelet_type,mode='symmetric', level=9, axis=-1)
    
    DWTcoeffs[-1] = np.zeros_like(DWTcoeffs[-1])
    DWTcoeffs[-2] = np.zeros_like(DWTcoeffs[-2])
    DWTcoeffs[-3] = np.zeros_like(DWTcoeffs[-3])
    DWTcoeffs[-4] = np.zeros_like(DWTcoeffs[-4])
    DWTcoeffs[-5] = np.zeros_like(DWTcoeffs[-5])
    DWTcoeffs[-6] = np.zeros_like(DWTcoeffs[-6])
    DWTcoeffs[-7] = np.zeros_like(DWTcoeffs[-7])
    DWTcoeffs[-8] = np.zeros_like(DWTcoeffs[-8])
    DWTcoeffs[-9] = np.zeros_like(DWTcoeffs[-9])
    
    filtered_data_dwt=pywt.waverec(DWTcoeffs,wavelet_type,mode='symmetric',axis=-1) 
    
    
    fc = 0.35  # Cut-off frequency of the butterworth filter
    w = fc / (360 / 2) # Normalize the frequency
    b, a = signal.butter(5, w, 'low')
    filtered_data_butterworth = signal.filtfilt(b, a, data)
    

    Let's plot the power-spectral density for the original data and the filtered data using both approaches:

    plt.figure(1)
    plt.psd(data, NFFT=512, Fs=360, label='original data', color='blue')
    plt.psd(filtered_data_dwt, NFFT=512, Fs=360, color='red', label='filtered data (DWT)')
    plt.psd(filtered_data_butterworth, NFFT=512, Fs=360, color='black', label='filtered data (Butterworth)')
    plt.legend()
    

    Which yields:

    enter image description here

    In the original data you can clearly see the 60Hz and its first multiple (120Hz). Let's take a close-up view at the low frequencies:

    enter image description here

    Now let's take a look at the data in time domain:

    plt.figure(2)
    plt.subplot(311)
    plt.plot(data,label='original data', color='blue')
    plt.title('original')
    plt.subplot(312)
    plt.plot(filtered_data_dwt, color='red', label='filtered data (DWT)')
    plt.title('filtered (DWT)')
    plt.subplot(313)
    plt.plot(filtered_data_butterworth, color='black', label='filtered data (Butterworth)')
    plt.title('filtered (Butterworth)')
    

    enter image description here

    So in order to low-pass filter your original data using a cut-off frequency of 0.35Hz, you can simply reconstruct them using the approximation coefficients of the DWT decomposition (that is, using a 'db6' wavelet). Hope this helps!