Как совместить вейвлет-преобразование и частотную фильтрацию

Question

Мне нужно реализовать следующее подавление шума на ЭКГ-сигнале:

• Дискретное вейвлет-преобразование до 9 уровней с помощью вейвлета 'db6'
• Фильтровать частоты (не детализировать коэффициенты)) на 9-м уровне в диапазоне 0–0,35 Гц
• Реконструировать сигнал, используя только уровни от 3 до 9

Я не знаю, как выполнить второй шаг вPython (PyWavelets), потому что я могу изменять только детали и коэффициенты аппроксимации, и я не знаю, как соотнести их с частотами.

Как мне поступить?

Это мой код

    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' )

Answer 1

Мой предыдущий (теперь удаленный) ответ был немного запутанным.Здесь я постараюсь предоставить вам практический пример, показывающий, что восстановление данных ЭКГ, отобранных на частоте 360 Гц с использованием только коэффициентов аппроксимации db6, (приблизительно) эквивалентно фильтрации нижних частот этих данных с использованием частоты среза 0,35 Гц.

В приведенном ниже примере кода я импортирую временной ряд ЭКГ из scipy (from scipy.misc import electrocardiogram);они сэмплированы на частоте 360 Гц, как и у вас.Я отфильтрую эти данные, используя:

• подход DWT, то есть реконструкцию данных с использованием только коэффициентов аппроксимации (Filter_Data_DWT)
• фильтра Баттерворта (Filter_Data_Butterworth)

Вот пример кода:

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)

Построим спектральную плотность мощности для исходных данных и отфильтрованных данных, используя оба подхода:

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()

Что дает:

В исходных данных четко видны 60 Гц и его первое кратное (120 Гц).Давайте рассмотрим крупным планом на низких частотах:

Теперь давайте посмотрим на данные во временной области:

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)')

Таким образом, чтобы выполнить низкочастотную фильтрацию исходных данных, используя частоту среза 0,35 Гц, вы можете просто восстановить их, используяаппроксимационные коэффициенты разложения DWT (то есть с использованием вейвлета 'db6'). Надеюсь, это поможет!