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

问题描述

我正在尝试对从采样频率为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应用了最低平滑度，但这也无法提供令人满意的结果.

解决方案

以下是一些建议.

首先，使用it=0尝试statsmodels中的lowess函数，并稍微调整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>]

第二个建议是使用 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()

这是结果的图.注意右边缘滤波信号的偏差.为了解决这个问题，您可以尝试使用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屋！