wake-up-neo.net

Wie kann ich mit SciPy/Numpy filtern/glätten?

Ich versuche, das von einem Druckwandler mit einer Abtastfrequenz von 50 kHz erhaltene Signal zu filtern/zu glätten. Ein Beispielsignal ist unten gezeigt:

enter image description here

Ich möchte ein glattes Signal erhalten, das durch Löss in MATLAB erhalten wird (ich zeichne nicht die gleichen Daten auf, die Werte sind unterschiedlich).

enter image description here

Ich habe die spektrale Leistungsdichte mit der psd () - Funktion von matplotlib berechnet, und die spektrale Leistungsdichte ist unten angegeben:

enter image description here

Ich habe versucht, folgenden Code zu verwenden und ein gefiltertes Signal zu erhalten:

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

Die Ausgabe, die ich bekomme, ist:

enter image description here

Ich brauche mehr Glättung, ich habe versucht, die Grenzfrequenz zu ändern, aber immer noch keine zufriedenstellenden Ergebnisse erhalten. Ich kann mit MATLAB nicht die gleiche Glätte bekommen. Ich bin sicher, dass dies in Python möglich ist, aber wie?

Sie finden die Daten hier .

Update

Ich habe eine geringe Glättung von Statistikmodellen angewendet, dies liefert ebenfalls keine zufriedenstellenden Ergebnisse.

enter image description here

pythonnumpyscipyfilteringsmoothing
15
Nimal Naser

Hier sind ein paar Vorschläge.

Probieren Sie zunächst die lowess-Funktion von statsmodels mit it=0 aus und verändern Sie das frac-Argument etwas:

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

Ein zweiter Vorschlag ist die Verwendung von scipy.signal.filtfilt anstelle von lfilter, um den Butterworth-Filter anzuwenden. filtfilt ist der Filter forward-backward. Es wendet den Filter zweimal an, einmal vorwärts und einmal rückwärts, was zu einer Phasenverzögerung von Null führt.

Hier ist eine modifizierte Version Ihres Skripts. Die wesentlichen Änderungen sind die Verwendung von filtfilt anstelle von lfilter und die Änderung von cutoff von 3000 auf 1500. Möglicherweise möchten Sie mit diesem Parameter experimentieren. Höhere Werte führen zu einer besseren Verfolgung des Beginns des Druckanstiegs, jedoch zu einem Wert zu hoch filtert die 3kHz (ungefähr) Schwingungen nach dem Druckanstieg nicht aus.

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

Hier ist die Darstellung des Ergebnisses. Beachten Sie die Abweichung des gefilterten Signals am rechten Rand. Um dies zu erreichen, können Sie mit den Parametern padtype und padlen von filtfilt experimentieren. Wenn Sie wissen, dass Sie über genügend Daten verfügen, können Sie die Flanken des gefilterten Signals verwerfen.

plot of filtfilt result

23
Warren Weckesser

Hier ist ein Beispiel für die Verwendung einer loewess-Anpassung. Beachten Sie, dass ich den Header von data.dat entfernt habe.

Aus der Darstellung scheint es, dass diese Methode bei Teilmengen der Daten gut funktioniert. Das Zusammenfügen aller Daten ergibt kein vernünftiges Ergebnis. Also ist dies wahrscheinlich nicht die beste Methode.

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

Das ist das Ergebnis, das ich bekomme:

plot

0
cel