如何使用 SciPy/Numpy 过滤/平滑? [英] How to filter/smooth with SciPy/Numpy?

问题描述

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

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

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

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

导入csv将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt将 scipy 导入为 sp从 scipy.signal 导入黄油、lfilter、freqzdef butter_lowpass(cutoff, fs, order=5):nyq = 0.5 * fsnormal_cutoff = 截止/nyqb, a = butter(order, normal_cutoff, btype='low',analog=False)返回 b, adef butter_lowpass_filter(data, cutoff, fs, order=5):b, a = butter_lowpass(cutoff, fs, order=order)y = lfilter(b, a, 数据)返回 ydata = np.loadtxt('data.dat', skiprows=2, delimiter=',', unpack=True).transpose()时间 = 数据[:,0]压力 = 数据[:,1]截止值 = 2000fs = 50000pressure_smooth = butter_lowpass_filter(压力，截止，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(时间，压力，线宽=1.0)plot_P_vs_t.plot(时间，压力平滑，线宽=1.0)plot_P_vs_t.set_ylabel('压力(bar)', labelpad=6)plot_P_vs_t.set_xlabel('时间(毫秒)'，labelpad=6)plt.show()plt.close()

我得到的输出是:

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

您可以在 而不是 lfilter 来应用巴特沃斯过滤器.filtfilt 是 forward-backward 过滤器.它应用了两次滤波器，一次向前，一次向后，导致零相位延迟.

这是您脚本的修改版本.显着的变化是使用 filtfilt 而不是 lfilter，以及将 cutoff 从 3000 更改为 1500.您可能想尝试一下参数 -- 较高的值会更好地跟踪压力增加的开始，但过高的值不会过滤掉压力增加后的 3kHz(大致)振荡.

将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt从 scipy.signal 进口黄油，filtfiltdef butter_lowpass(cutoff, fs, order=5):nyq = 0.5 * fsnormal_cutoff = 截止/nyqb, a = butter(order, normal_cutoff, btype='low',analog=False)返回 b, adef butter_lowpass_filtfilt(data, cutoff, fs, order=5):b, a = butter_lowpass(cutoff, fs, order=order)y = filtfilt(b, a, 数据)返回 ydata = np.loadtxt('data.dat', skiprows=2, delimiter=',', unpack=True).transpose()时间 = 数据[:,0]压力 = 数据[:,1]截止值 = 1500fs = 50000pressure_smooth = butter_lowpass_filtfilt(压力，截止，fs)figure_pressure_trace = plt.figure()figure_pressure_trace.clf()plot_P_vs_t = plt.subplot(111)plot_P_vs_t.plot(时间，压力，'r'，线宽=1.0)plot_P_vs_t.plot(时间，压力平滑，'b'，线宽=1.0)plt.show()plt.close()

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

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

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

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

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?

You can find the data here.

Update

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

解决方案

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.

这篇关于如何使用 SciPy/Numpy 过滤/平滑?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！