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

Question

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

我想在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应用了lowess平滑,这也没有提供令人满意的结果.

Answer 1

以下是一些建议.

首先,尝试使用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>]

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

这是您脚本的修改版本.显着的变化是使用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()

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