代码之家 › 专栏 › 技术社区 › Nimal Naser

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

smoothing filtering scipy numpy python

Nimal Naser · 技术社区 · 11 年前

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

enter image description here

我希望在MATLAB中获得黄土获得的平滑信号(我绘制的不是相同的数据,值不同)。

enter image description here

我使用matplotlib的psd()函数计算了功率谱密度,功率谱密度也如下所示:

enter image description here

我已尝试使用以下代码,并获得了一个滤波信号:

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

我得到的输出是:

enter image description here

我需要更多的平滑,我尝试改变截止频率,但仍然无法获得满意的结果。我无法通过MATLAB获得相同的平滑度。我相信这可以在Python中实现,但如何实现?

您可以找到数据 here .

使现代化

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

enter image description here

2 回复 | 直到 11 年前

Warren Weckesser 11 年前

这里有几个建议。

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

plot

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

这是您的脚本的修改版本。重要的变化是使用 过滤,过滤 而不是 硫过滤器 ,以及 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 的参数 过滤,过滤 或者,如果你知道你有足够的数据,你可以丢弃滤波信号的边缘。

plot of filtfilt result

cel 11 年前

下面是使用 loewess 适合注意,我从 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()

这是我得到的结果:

plot