python - functions - Cómo filtrar/suavizar con SciPy/Numpy?

numpy windows

Question

Estoy tratando de filtrar / suavizar la señal obtenida de un transductor de presión de frecuencia de muestreo de 50 kHz. Una señal de muestra se muestra a continuación:

Me gustaría obtener una señal suave obtenida por loess en MATLAB (no estoy trazando los mismos datos, los valores son diferentes).

Calculé la densidad espectral de potencia usando la función psd () de matplotlib y la densidad espectral de potencia también se proporciona a continuación:

Intenté usar el siguiente código y obtuve una señal filtrada:

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

El resultado que obtengo es:

Necesito más suavizado, intenté cambiar la frecuencia de corte, pero todavía no se pueden obtener resultados satisfactorios. No puedo obtener la misma suavidad de MATLAB. Estoy seguro de que se puede hacer en Python, pero ¿cómo?

Puede encontrar los datos aquí .

Actualizar

Apliqué suavizar lowess desde modelos de estadísticas, esto tampoco proporciona resultados satisfactorios.

Answer 1

Aquí hay un ejemplo de usar un ajuste loewess . Tenga en cuenta que data.dat el encabezado de data.dat .

De la trama, parece que este método funciona bien en subconjuntos de datos. Ajustar todos los datos a la vez no da un resultado razonable. Entonces, probablemente este no sea el mejor método aquí.

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

Este es el resultado que obtengo:

Answer 2

Aquí hay un par de sugerencias.

Primero, pruebe la función lowess de los statsmodels de statsmodels con it=0 , y statsmodels un poco el argumento 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>]

Una segunda sugerencia es usar scipy.signal.filtfilt lugar de lfilter para aplicar el filtro Butterworth. filtfilt es el filtro hacia adelante-hacia atrás . Aplica el filtro dos veces, una vez hacia adelante y otra hacia atrás, lo que produce un retraso de fase cero.

Aquí hay una versión modificada de su script. Los cambios significativos son el uso de filtfilt lugar de lfilter y el cambio de cutoff de 3000 a 1500. Es posible que desee experimentar con ese parámetro: valores más altos dan como resultado un mejor seguimiento del inicio del aumento de presión, pero un valor que es demasiado alto no filtra las oscilaciones de 3kHz (aproximadamente) después de que la presión aumenta.

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

Aquí está la trama del resultado. Tenga en cuenta la desviación en la señal filtrada en el borde derecho. Para manejar eso, puede experimentar con los parámetros padtype y padlen de filtfilt o, si sabe que tiene suficientes datos, puede descartar los bordes de la señal filtrada.