Как фильтровать/сглаживать SciPy/Numpy?

Я пытаюсь фильтровать/сглаживать сигнал, полученный от датчика давления частоты дискретизации 50 кГц. Пример сигнала показан ниже:

Я хотел бы получить плавный сигнал, полученный лёссой в MATLAB (я не рисую одни и те же данные, значения разные).

Я вычислил спектральную плотность мощности, используя функцию matplotlib psd(), и спектральная плотность мощности также представлена ниже:

Я попытался использовать следующий код и получил отфильтрованный сигнал:

import csv
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.signal import butter, lfilter, freqz

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

data = np.loadtxt('data.dat', skiprows=2, delimiter=',', unpack=True).transpose()
time = data[:,0]
pressure = data[:,1]
cutoff = 2000
fs = 50000
pressure_smooth = butter_lowpass_filter(pressure, cutoff, fs)

figure_pressure_trace = plt.figure(figsize=(5.15, 5.15))
figure_pressure_trace.clf()
plot_P_vs_t = plt.subplot(111)
plot_P_vs_t.plot(time, pressure, linewidth=1.0)
plot_P_vs_t.plot(time, pressure_smooth, linewidth=1.0)
plot_P_vs_t.set_ylabel('Pressure (bar)', labelpad=6)
plot_P_vs_t.set_xlabel('Time (ms)', labelpad=6)
plt.show()
plt.close()

Выход, который я получаю:

Мне нужно больше сглаживания, я попытался изменить частоту среза, но все же удовлетворительные результаты не могут быть получены. MATLAB не может получить такую же гладкость. Я уверен, что это можно сделать на Python, но как?

Здесь вы можете найти данные .

Обновление

Я применил понижающее сглаживание от statsmodels, это также не дает удовлетворительных результатов.

Ответ 1

Вот пара предложений.

Сначала попробуйте функцию lowess от statsmodels с помощью it=0 и немного измените аргумент frac:

In [328]: from statsmodels.nonparametric.smoothers_lowess import lowess

In [329]: filtered = lowess(pressure, time, is_sorted=True, frac=0.025, it=0)

In [330]: plot(time, pressure, 'r')
Out[330]: [<matplotlib.lines.Line2D at 0x1178d0668>]

In [331]: plot(filtered[:,0], filtered[:,1], 'b')
Out[331]: [<matplotlib.lines.Line2D at 0x1173d4550>]

Второе предложение - использовать scipy.signal.filtfilt вместо lfilter для применения фильтра Баттерворта. filtfilt - фильтр с обратной связью. Он применяет фильтр дважды, один раз вперед и один раз назад, что приводит к задержке нулевой фазы.

Здесь представлена измененная версия вашего script. Существенными изменениями являются использование filtfilt вместо lfilter и изменение cutoff от 3000 до 1500. Возможно, вам захочется поэкспериментировать с этим параметром - более высокие значения приводят к лучшему отслеживанию начала давления но слишком высокое значение не отфильтровывает 3 кГц (грубые) колебания после повышения давления.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

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_filtfilt(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)
    return y

data = np.loadtxt('data.dat', skiprows=2, delimiter=',', unpack=True).transpose()
time = data[:,0]
pressure = data[:,1]
cutoff = 1500
fs = 50000
pressure_smooth = butter_lowpass_filtfilt(pressure, cutoff, fs)

figure_pressure_trace = plt.figure()
figure_pressure_trace.clf()
plot_P_vs_t = plt.subplot(111)
plot_P_vs_t.plot(time, pressure, 'r', linewidth=1.0)
plot_P_vs_t.plot(time, pressure_smooth, 'b', linewidth=1.0)
plt.show()
plt.close()

Вот график результата. Обратите внимание на отклонение отфильтрованного сигнала на правом краю. Чтобы справиться с этим, вы можете поэкспериментировать с параметрами padtype и padlen filtfilt, или, если вы знаете, что у вас достаточно данных, вы можете отбросить края отфильтрованного сигнала.

Ответ 2

Вот пример использования loewess fit. Обратите внимание, что я разделил заголовок с data.dat.

Из графика видно, что этот метод хорошо работает на подмножествах данных. Приведение всех данных одновременно не дает разумного результата. Итак, вероятно, это не лучший метод здесь.

import pandas as pd
import matplotlib.pylab as plt
from statsmodels.nonparametric.smoothers_lowess import lowess

data = pd.read_table("data.dat", sep=",", names=["time", "pressure"])
sub_data = data[data.time > 21.5]

result = lowess(sub_data.pressure, sub_data.time.values)
x_smooth = result[:,0]
y_smooth = result[:,1]

tot_result = lowess(data.pressure, data.time.values, frac=0.1)
x_tot_smooth = tot_result[:,0]
y_tot_smooth = tot_result[:,1]

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(data.time.values, data.pressure, label="raw")
ax.plot(x_tot_smooth, y_tot_smooth, label="lowess 1%", linewidth=3, color="g")
ax.plot(x_smooth, y_smooth, label="lowess", linewidth=3, color="r")
plt.legend()

В результате я получаю: