Làm thế nào để lọc/trơn tru với SciPy/Numpy?

Question

Tôi đang cố gắng lọc/tín hiệu trơn tru thu được từ bộ chuyển đổi áp suất của tần số lấy mẫu 50 kHz. Một tín hiệu mẫu được hiển thị dưới đây:

Tôi muốn để có được một tín hiệu mịn thu được bằng cách hoàng thổ trong MATLAB (Tôi không âm mưu cùng một dữ liệu, giá trị này là khác nhau).

tôi tính toán mật độ phổ công suất sử dụng chức năng matplotlib của psd() và mật độ phổ công suất cũng được cung cấp dưới đây:

Tôi đã thử bằng cách sử dụng đoạn mã sau và thu được một tín hiệu được lọc:

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

Sản lượng tôi nhận được là:

tôi cần làm mịn hơn, tôi cố gắng thay đổi tần số cắt nhưng vẫn kết quả khả quan không thể có được. Tôi không thể có được sự mượt mà như nhau bởi MATLAB. Tôi chắc chắn nó có thể được thực hiện bằng Python, nhưng làm thế nào?

Bạn có thể tìm thấy dữ liệu here.

Cập nhật

tôi áp dụng lowess mịn từ statsmodels, điều này cũng không cung cấp kết quả khả quan.

Answer 1

Dưới đây là một vài gợi ý.

Đầu tiên, hãy thử các lowess chức năng từ statsmodels với it=0, và tinh chỉnh các frac tranh luận một chút:

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>] 

Một gợi ý thứ hai là sử dụng scipy.signal.filtfilt thay vì lfilter áp dụng bộ lọc Butterworth . filtfilt là bộ lọc chuyển tiếp lùi lại. Nó áp dụng các bộ lọc hai lần, một lần về phía trước và một lần lạc hậu, dẫn đến sự chậm trễ pha số không.

Dưới đây là một phiên bản sửa đổi của kịch bản của bạn. Những thay đổi đáng kể là việc sử dụng các filtfilt thay vì lfilter, và sự thay đổi của cutoff từ 3000 đến 1500. Bạn có thể muốn thử nghiệm với tham số đó - giá trị cao hơn kết quả trong việc theo dõi tốt hơn về sự khởi đầu của sự gia tăng áp lực, nhưng là một giá trị mà quá cao không lọc ra dao động 3kHz (gần) sau khi tăng áp suất.

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

Đây là cốt truyện của kết quả. Lưu ý độ lệch trong tín hiệu được lọc ở cạnh phải. Để xử lý, bạn có thể thử nghiệm với các padtype và padlen thông số của filtfilt, hoặc, nếu bạn biết bạn có đủ dữ liệu, bạn có thể loại bỏ các cạnh của tín hiệu được lọc.

Answer 2

Dưới đây là một ví dụ của việc sử dụng một sự phù hợp loewess. Lưu ý rằng tôi đã tước tiêu đề từ data.dat.

Từ cốt truyện có vẻ như phương pháp này hoạt động tốt trên các tập con của dữ liệu. Việc lắp tất cả dữ liệu cùng một lúc không mang lại kết quả hợp lý. Vì vậy, có lẽ đây không phải là phương pháp tốt nhất ở đây.

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

Đây là kết quả tôi nhận được: