Python 如何使用 SciPy/Numpy 过滤/平滑？

How to filter/smooth with SciPy/Numpy?

提问by Nimal Naser

I am trying to filter/smooth signal obtained from a pressure transducer of sampling frequency 50 kHz. A sample signal is shown below:

我正在尝试过滤/平滑从采样频率为 50 kHz 的压力传感器获得的信号。示例信号如下所示：

I would like to obtain a smooth signal obtained by loess in MATLAB (I am not plotting the same data, values are different).

我想在 MATLAB 中获得由 loess 获得的平滑信号（我没有绘制相同的数据，值不同）。

I calculated the power spectral density using matplotlib's psd() function and the power spectral density is also provided below:

我使用 matplotlib 的 psd() 函数计算了功率谱密度，下面还提供了功率谱密度：

I have tried using the following code and obtained a filtered signal:

我尝试使用以下代码并获得过滤信号：

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

The output I get is:

我得到的输出是：

I need more smoothing, I tried changing the cutoff frequency but still satisfactory results can not be obtained. I can't get the same smoothness by MATLAB. I am sure it can be done in Python, but how?

我需要更多的平滑，我尝试改变截止频率，但仍然无法获得令人满意的结果。我无法通过 MATLAB 获得相同的平滑度。我确信它可以在 Python 中完成，但是如何呢？

You can find the data here.

您可以在此处找到数据。

Update

更新

I applied lowess smoothing from statsmodels, this also does not provide satisfactory results.

我从 statsmodels 应用了 Lowess 平滑，这也不能提供令人满意的结果。

采纳答案by Warren Weckesser

Here are a couple suggestions.

这里有一些建议。

First, try the lowessfunction from statsmodelswith it=0, and tweak the fracargument a bit:

首先，尝试使用with的lowess函数，并稍微调整参数：statsmodelsit=0frac

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.filtfiltinstead of lfilterto apply the Butterworth filter. filtfiltis the forward-backwardfilter. It applies the filter twice, once forward and once backward, resulting in zero phase delay.

第二个建议是使用scipy.signal.filtfilt而不是lfilter应用巴特沃斯滤波器。 filtfilt是前向后向滤波器。它应用滤波器两次，一次向前，一次向后，导致零相位延迟。

Here's a modified version of your script. The significant changes are the use of filtfiltinstead of lfilter, and the change of cutofffrom 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.

这是您的脚本的修改版本。显着的变化是使用filtfilt代替lfilter，以及cutoff从 3000 到 1500的变化。您可能想用该参数进行试验——较高的值会导致更好地跟踪压力增加的开始，但值太高在压力增加后不会过滤掉 3kHz（大致）振荡。

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 padtypeand padlenparameters of filtfilt, or, if you know you have enough data, you can discard the edges of the filtered signal.

这是结果图。注意右边缘滤波信号的偏差。为了解决这个问题，您可以试验 的padtype和padlen参数filtfilt，或者，如果您知道自己有足够的数据，则可以丢弃滤波信号的边缘。

回答by cel

Here is an example of using a loewessfit. Note that I stripped the header from data.dat.

这是使用loewess拟合的示例。请注意，我从data.dat.

From the plot it seems that this method performs well on subsets of the data. Fitting all data at once does not give a reasonable result. So, probably this is not the best method here.

从图中看来，这种方法在数据子集上表现良好。一次拟合所有数据并不能给出合理的结果。所以，这可能不是最好的方法。

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

This is the result I get:

这是我得到的结果：