¿Cómo filtrar / suavizar con SciPy / Numpy?

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

introduzca la descripción de la imagen aquí

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

introduzca la descripción de la imagen aquí

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

    introduzca la descripción de la imagen aquí

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

    La salida que obtengo es:

    introduzca la descripción de la imagen aquí

    Necesito más suavizado, intenté cambiar la frecuencia de corte pero aún 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?

    Puedes encontrar los datos aquí .

    Actualizar

    Apliqué lowess smoothing desde statsmodels, esto tampoco proporciona resultados satisfactorios.

    introduzca la descripción de la imagen aquí

    Aquí hay un par de sugerencias.

    Primero, pruebe la función lowess de statsmodels con it=0 , y statsmodels un poco el argumento de la statsmodels :

     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]: [] In [331]: plot(filtered[:,0], filtered[:,1], 'b') Out[331]: [] 

    trama

    Una segunda sugerencia es usar scipy.signal.filtfilt lugar de lfilter para aplicar el filtro Butterworth. filtfilt es el filtro de avance hacia atrás . Aplica el filtro dos veces, una vez hacia adelante y una 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: los valores más altos resultan en un mejor seguimiento del inicio del aumento de presión, pero un valor que es demasiado alto no filtra las oscilaciones de 3 kHz (aproximadamente) después del aumento de presión.

     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 ttwig del resultado. Note 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.

    trama del resultado de filtfilt

    Aquí hay un ejemplo del uso de un ajuste de loewess . Tenga en cuenta que he quitado el encabezado de data.dat .

    De la ttwig parece que este método funciona bien en subconjuntos de los 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:

    trama