且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

如何使用 SciPy/Numpy 过滤/平滑?

更新时间:2023-02-26 16:41:55

Here are a couple suggestions.

First, try the lowess function from statsmodels with it=0, and tweak the frac argument a bit:

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>]

A second suggestion is to use scipy.signal.filtfilt instead of lfilter to apply the Butterworth filter. filtfilt is the forward-backward filter. It applies the filter twice, once forward and once backward, resulting in zero phase delay.

Here's a modified version of your script. The significant changes are the use of filtfilt instead of lfilter, and the change of cutoff from 3000 to 1500. You might want to experiment with that parameter--higher values result in better tracking of the onset of the pressure increase, but a value that is too high doesn't filter out the 3kHz (roughly) oscillations after the pressure increase.

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

Here's the plot of the result. Note the deviation in the filtered signal at the right edge. To handle that, you can experiment with the padtype and padlen parameters of filtfilt, or, if you know you have enough data, you can discard the edges of the filtered signal.