¿Cómo suavizar una curva de la manera correcta?

200

Supongamos que tenemos un conjunto de datos que podría ser dado aproximadamente por

import numpy as np
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2

Por lo tanto, tenemos una variación del 20% del conjunto de datos. Mi primera idea fue usar la función UnivariateSpline de scipy, pero el problema es que esto no considera el ruido pequeño en el buen sentido. Si considera las frecuencias, el fondo es mucho más pequeño que la señal, por lo que una spline solo del corte podría ser una idea, pero eso implicaría una transformación de Fourier de ida y vuelta, lo que podría provocar un mal comportamiento. Otra forma sería un promedio móvil, pero esto también necesitaría la elección correcta del retraso.

¿Alguna pista / libro o enlace sobre cómo abordar este problema?

ejemplo

varantir
fuente
1
¿Su señal siempre será una onda sinusoidal, o la estaba usando solo como ejemplo?
Mark Ransom
no, tendré diferentes señales, incluso en este sencillo ejemplo, es obvio que mis métodos no son suficientes
varantir
El filtrado Kalman es óptimo para este caso. Y el paquete pykalman python es de buena calidad.
toine
Tal vez lo amplíe a una respuesta completa cuando tenga un poco más de tiempo, pero el único método de regresión poderoso que aún no se mencionó es la regresión GP (Proceso Gaussiano).
Ori5678

Respuestas:

262

Prefiero un filtro Savitzky-Golay . Utiliza mínimos cuadrados para hacer retroceder una pequeña ventana de sus datos en un polinomio, luego usa el polinomio para estimar el punto en el centro de la ventana. Finalmente, la ventana se desplaza hacia adelante por un punto de datos y el proceso se repite. Esto continúa hasta que cada punto se haya ajustado de manera óptima en relación con sus vecinos. Funciona muy bien incluso con muestras ruidosas de fuentes no periódicas y no lineales.

Aquí hay un ejemplo completo de libro de cocina . Vea mi código a continuación para tener una idea de lo fácil que es usarlo. Nota: omití el código para definir la savitzky_golay()función porque literalmente puedes copiarlo / pegarlo del ejemplo del libro de cocina que vinculé anteriormente.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
yhat = savitzky_golay(y, 51, 3) # window size 51, polynomial order 3

plt.plot(x,y)
plt.plot(x,yhat, color='red')
plt.show()

alisando óptimamente una sinusoide ruidosa

ACTUALIZACIÓN: Me ha llamado la atención que el ejemplo del libro de cocina al que me vinculé ha sido eliminado. Afortunadamente, el filtro Savitzky-Golay se ha incorporado a la biblioteca SciPy , como lo señaló @dodohjk . Para adaptar el código anterior utilizando la fuente SciPy, escriba:

from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
David Wurtz
fuente
Recibí el error Traceback (última llamada): archivo "hp.py", línea 79, en <module> ysm2 = savitzky_golay (y_data, 51,3) Archivo "hp.py", línea 42, en savitzky_golay firstvals = y [0] - np.abs (y [1: half_window + 1] [:: - 1] - y [0])
Marzo Ho
14
¡Gracias por presentar el filtro Savitzky-Golay! Básicamente, esto es como un filtro normal de "promedio móvil", pero en lugar de calcular el promedio, se hace un ajuste polinómico (generalmente de segundo o cuarto orden) para cada punto, y solo se elige el punto "medio". Dado que la información de segundo (o cuarto) orden se refiere a cada punto, se evita el sesgo introducido en el enfoque de "promedio móvil" en máximos o mínimos locales. Muy elegante
np8
2
Solo quiero decir gracias por esto, me he vuelto loco tratando de descubrir las descomposiciones de wavelet para obtener datos suavizados, y esto es mucho mejor.
Eldar M.
55
Si los datos x no está regularmente espaciados es posible que desee aplicar el filtro de las x, así: savgol_filter((x, y), ...).
Tim Kuipers
127

Una forma rápida y sucia de suavizar los datos que uso, basada en un cuadro de media móvil (por convolución):

x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.8

def smooth(y, box_pts):
    box = np.ones(box_pts)/box_pts
    y_smooth = np.convolve(y, box, mode='same')
    return y_smooth

plot(x, y,'o')
plot(x, smooth(y,3), 'r-', lw=2)
plot(x, smooth(y,19), 'g-', lw=2)

ingrese la descripción de la imagen aquí

scrx2
fuente
9
Esto tiene algunas ventajas interesantes: (1) funciona para cualquier función, no solo periódica, y (2) no depende de funciones grandes para copiar y pegar. Puedes hacerlo de inmediato con Numpy puro. Además, no está demasiado sucio: es el caso más simple de algunos de los otros métodos descritos anteriormente (como BAJO pero el núcleo es un intervalo agudo y como Savitzky-Golay pero el grado polinómico es cero).
Jim Pivarski
2
El único problema con la media móvil es que se queda atrás de los datos. Puede ver esto más aparentemente al final, donde hay más puntos en la parte superior y menos en la parte inferior, pero la curva verde está actualmente por debajo del promedio porque la función de ventana tiene que avanzar para tenerlos en cuenta.
nurettin
Y esto no funciona en nd array, solo 1d. scipy.ndimage.filters.convolve1d()le permite especificar un eje de una matriz nd para hacer el filtrado. Pero creo que ambos sufren algunos problemas en los valores enmascarados.
Jason
1
@nurettin Creo que lo que estás describiendo son efectos de borde. En general, siempre que el núcleo de convolución pueda cubrir su extensión dentro de la señal, no se "queda atrás" como usted dice. Al final, sin embargo, no hay valores más allá de 6 para incluir en el promedio, por lo que solo se está utilizando la porción "izquierda" del núcleo. Los efectos de borde están presentes en cada núcleo de suavizado y deben manejarse por separado.
Jon
44
@nurettin No, estaba tratando de aclararle a los demás que leen esto que su comentario "el único problema con la media móvil es que va a la zaga de los datos" es engañoso. Cualquier método de filtro de ventana sufre este problema, no solo el promedio móvil. Savitzky-golay también sufre este problema. Entonces su afirmación "Lo que estoy describiendo es lo que savitzky_golay resuelve por estimación" es simplemente errónea. Cualquiera de los métodos de suavizado requiere una forma de manejar bordes que sea independiente del método de suavizado en sí.
Jon
79

Si está interesado en una versión "uniforme" de una señal que sea periódica (como su ejemplo), entonces una FFT es el camino correcto. Tome la transformada de Fourier y reste las frecuencias de baja contribución:

import numpy as np
import scipy.fftpack

N = 100
x = np.linspace(0,2*np.pi,N)
y = np.sin(x) + np.random.random(N) * 0.2

w = scipy.fftpack.rfft(y)
f = scipy.fftpack.rfftfreq(N, x[1]-x[0])
spectrum = w**2

cutoff_idx = spectrum < (spectrum.max()/5)
w2 = w.copy()
w2[cutoff_idx] = 0

y2 = scipy.fftpack.irfft(w2)

ingrese la descripción de la imagen aquí

Incluso si su señal no es completamente periódica, esto hará un gran trabajo restando el ruido blanco. Hay muchos tipos de filtros para usar (paso alto, paso bajo, etc.), el apropiado depende de lo que esté buscando.

Enganchado
fuente
¿Qué parcela es para qué variable? Estoy tratando de suavizar las coordenadas de la pelota de tenis en un rally, es decir. sacar todos los rebotes que parecen pequeñas parábolas en mi trama
mLstudent33
44

Ajustar un promedio móvil a sus datos suavizaría el ruido, vea esta respuesta para saber cómo hacerlo.

Si desea utilizar LOWESS para ajustar sus datos (es similar a un promedio móvil pero más sofisticado), puede hacerlo utilizando la biblioteca de statsmodels :

import numpy as np
import pylab as plt
import statsmodels.api as sm

x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
lowess = sm.nonparametric.lowess(y, x, frac=0.1)

plt.plot(x, y, '+')
plt.plot(lowess[:, 0], lowess[:, 1])
plt.show()

Finalmente, si conoce la forma funcional de su señal, podría ajustar una curva a sus datos, lo que probablemente sería lo mejor.

markmuetz
fuente
Si solo la hubiera loessimplementado.
scrutari
18

Otra opción es usar KernelReg en statsmodels :

from statsmodels.nonparametric.kernel_regression import KernelReg
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2

# The third parameter specifies the type of the variable x;
# 'c' stands for continuous
kr = KernelReg(y,x,'c')
plt.plot(x, y, '+')
y_pred, y_std = kr.fit(x)

plt.plot(x, y_pred)
plt.show()
Zichen Wang
fuente
7

¡Mira esto! Hay una definición clara de suavizado de una señal 1D.

http://scipy-cookbook.readthedocs.io/items/SignalSmooth.html

Atajo:

import numpy

def smooth(x,window_len=11,window='hanning'):
    """smooth the data using a window with requested size.

    This method is based on the convolution of a scaled window with the signal.
    The signal is prepared by introducing reflected copies of the signal 
    (with the window size) in both ends so that transient parts are minimized
    in the begining and end part of the output signal.

    input:
        x: the input signal 
        window_len: the dimension of the smoothing window; should be an odd integer
        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
            flat window will produce a moving average smoothing.

    output:
        the smoothed signal

    example:

    t=linspace(-2,2,0.1)
    x=sin(t)+randn(len(t))*0.1
    y=smooth(x)

    see also: 

    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
    scipy.signal.lfilter

    TODO: the window parameter could be the window itself if an array instead of a string
    NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.
    """

    if x.ndim != 1:
        raise ValueError, "smooth only accepts 1 dimension arrays."

    if x.size < window_len:
        raise ValueError, "Input vector needs to be bigger than window size."


    if window_len<3:
        return x


    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"


    s=numpy.r_[x[window_len-1:0:-1],x,x[-2:-window_len-1:-1]]
    #print(len(s))
    if window == 'flat': #moving average
        w=numpy.ones(window_len,'d')
    else:
        w=eval('numpy.'+window+'(window_len)')

    y=numpy.convolve(w/w.sum(),s,mode='valid')
    return y




from numpy import *
from pylab import *

def smooth_demo():

    t=linspace(-4,4,100)
    x=sin(t)
    xn=x+randn(len(t))*0.1
    y=smooth(x)

    ws=31

    subplot(211)
    plot(ones(ws))

    windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']

    hold(True)
    for w in windows[1:]:
        eval('plot('+w+'(ws) )')

    axis([0,30,0,1.1])

    legend(windows)
    title("The smoothing windows")
    subplot(212)
    plot(x)
    plot(xn)
    for w in windows:
        plot(smooth(xn,10,w))
    l=['original signal', 'signal with noise']
    l.extend(windows)

    legend(l)
    title("Smoothing a noisy signal")
    show()


if __name__=='__main__':
    smooth_demo()
IPhysResearch
fuente
3
Un enlace a una solución es bienvenido, pero asegúrese de que su respuesta sea útil sin él: agregue contexto alrededor del enlace para que sus otros usuarios tengan una idea de qué es y por qué está allí, luego cite la parte más relevante de la página ' volver a vincular en caso de que la página de destino no esté disponible. Se pueden eliminar las respuestas que son poco más que un enlace.
Shree
-4

Si está trazando un gráfico de series de tiempo y si ha utilizado mtplotlib para dibujar gráficos, utilice el método mediano para suavizar el gráfico

smotDeriv = timeseries.rolling(window=20, min_periods=5, center=True).median()

¿Dónde se timeseriespasa su conjunto de datos que puede modificar windowsizepara obtener más suavizado?

Pavan Purohit
fuente