Применение подходящего фильтра Баттерворта к необработанному сигналу с использованием Python

Question

Я получил 10-секундный необработанный сигнал PPG (Фотоплетизмограмма) от моего TI AFE4490. Мое оборудование откалибровано, и я использую 250 выборок в секунду для записи этого сигнала. Я набрал 2500 очков в конце.

Я использовал полосовой фильтр Баттерворта с lowcut = 0.5, highcut = 15 и order = 2. Вы можете видеть мои необработанные и отфильтрованные сигналы ниже:

Я также пытался отфильтровать это, используя низкочастотный фильтр Баттерворта с lowcut = 15 и order = 2. Как видите, мои необработанные и отфильтрованные сигналы приведены ниже:

Я читал в некоторых статьях, что 0,5 Гц и 15 Гц являются хорошими низкими и высокими частотами для этого типа сигнала.

Перед тем, как применить фильтры, я использовал алгоритм Сципи Баттерворта (Scipy Docs), чтобы показать мне ответ фильтра, и это было хорошо.

Мой отфильтрованный сигнал кажется хорошим после этого «старта» (повышения в начале), но я не знаю, почему это начало. Кто-нибудь может сказать мне, если это "начало" это нормально на фильтрах Баттерворта? Если да, есть какой-нибудь способ исправить это?

Я ценю вашу помощь.

Мой код:

RED, IR, nSamples, sRate = getAFESignal()

period = 1/sRate # Signal period.
# Desired cutoff frequency (in Hz) and filter order.
lowcut = 0.5
highcut = 15
orders = 2

plt.figure(1)

x = np.linspace(0, nSamples*period, nSamples, endpoint=True)

plt.subplot(2,1,1)
y = IR
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.plot(x,y, label='Noisy signal')


plt.subplot(2,1,2)
yf = butter_bandpass_filter(IR, lowcut, highcut, nSamples, order=orders)
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.plot(x, yf, label='Filtered signal')

plt.grid()
plt.show()

Функция getAFEsignal() - это просто функция для чтения .txt файла и помещения всех в два массива.

Answer 1

Начальный переходный процесс, который вы видите на графиках, представляет собой шаговый отклик фильтра, когда внезапный вход применяется к фильтру в состоянии покоя. Если вы только что подключили физический прибор, включающий такой полосовой фильтр, датчик прибора мог бы получить выборки входных данных, прыгнув от 0 (пока зонд отключен) до значения первой выборки ~ 0,126 В. Отклик фильтра прибора показал бы аналогичный переходный процесс.

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

Один из способов добиться этого - использовать достаточно длинную выборку данных и отбрасывать начальный переходный процесс. Другой подход заключается в том, чтобы заставить начальное внутреннее состояние фильтра приблизиться к тому, что можно было бы ожидать, если бы сигнал аналогичной амплитуды подавался в течение некоторого времени перед вашей первой входной выборкой. Это можно сделать, например, установив начальное условие с помощью scipy.signal.lfilter_zi.

Теперь, я полагаю, вы использовали butter_bandpass_filter из SciPy Cookbook , который не заботится о начальных условиях фильтра. К счастью, это может быть легко изменено с этой целью:

from scipy.signal import butter, lfilter, lfilter_zi

def butter_bandpass_filter_zi(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    zi = lfilter_zi(b, a)
    y,zo = lfilter(b, a, data, zi=zi*data[0])
    return y

Еще один момент, на который следует обратить внимание, это то, что вы звоните butter_bandpass_filter как:

yf = butter_bandpass_filter(IR, lowcut, highcut, nSamples, order=orders)

Передача nSamples (общее количество выборок, в вашем случае 2500) в качестве 4-го аргумента, тогда как функция ожидает частоту дискретизации (в вашем случае 250). Коэффициент 10 между двумя величинами имеет эффект, эквивалентный уменьшению диапазона фильтрации с [0.5,15] Гц до [0.05,1.5] Гц. Чтобы получить предполагаемый диапазон частот полосы пропускания, вы должны передать sRate в качестве четвертого аргумента:

yf = butter_bandpass_filter_zi(IR, lowcut, highcut, sRate, order=orders)

Наконец, вы можете заметить, что этот последний вывод немного менее треугольный, чем вход. Это вызвано тем, что часть низкочастотного контента около 0,5 Гц отфильтровывается. Если это то, что вы ожидали, то отлично. В противном случае вы все равно можете поиграть с частотой среза, чтобы получить то, что, по вашему мнению, дает наилучшие результаты. Например (и я не хочу сказать, что это лучший частотный диапазон), если вы установите lowcut=0.25, вы получите более треугольный график, такой как: