Supongamos que tenemos un conjunto de datos que podría estar 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 pequeño ruido de una buena manera. 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 resultar en un mal comportamiento. Otra forma sería un promedio móvil, pero esto también necesitaría la elección correcta del retraso.

¿Algún consejo/libro o enlace sobre cómo abordar este problema?

prefiero un Filtro Savitzky-Golay. Utiliza mínimos cuadrados para 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 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 el savitzky_golay() función porque literalmente puede copiarlo/pegarlo del ejemplo del libro de cocina que vinculé arriba.

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

ACTUALIZAR: Me ha llamado la atención que el ejemplo del libro de cocina al que me vinculé ha sido eliminado. Afortunadamente se ha incorporado el filtro Savitzky-Golay en la biblioteca SciPy, como lo señaló @dodohjk (gracias @bicarlsen por el enlace actualizado). 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

EDITAR: mira esta respuesta. Usando np.cumsum es mucho más rápido que np.convolve

Una forma rápida y sucia de suavizar los datos que uso, basada en un cuadro de promedio 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)

Si está interesado en una versión “suave” de una señal que es 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)

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

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

Si desea utilizar BAJA para ajustar sus datos (es similar a un promedio móvil pero más sofisticado), puede hacerlo usando el modelos estadisticos biblioteca:

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 que podría hacer.

Esta pregunta ya está completamente respondida, por lo que creo que sería interesante un análisis en tiempo de ejecución de los métodos propuestos (lo fue para mí, de todos modos). También miraré el comportamiento de los métodos en el centro y los bordes del conjunto de datos ruidoso.

TL;DR

                    | runtime in s | runtime in s
method              | python list  | numpy array
--------------------|--------------|------------
kernel regression   | 23.93405     | 22.75967 
lowess              |  0.61351     |  0.61524 
naive average       |  0.02485     |  0.02326 
others*             |  0.00150     |  0.00150 
fft                 |  0.00021     |  0.00021 
numpy convolve      |  0.00017     |  0.00015 

*savgol with different fit functions and some numpy methods

La regresión del kernel escala mal, Lowess es un poco más rápido, pero ambos producen curvas suaves. Savgol es un término medio en cuanto a velocidad y puede producir salidas suaves y con saltos, según el grado del polinomio. FFT es extremadamente rápido, pero solo funciona con datos periódicos.

Los métodos de promedio móvil con numpy son más rápidos pero obviamente producen un gráfico con pasos.

Configuración

Generé 1000 puntos de datos en forma de curva sinusoidal:

size = 1000
x = np.linspace(0, 4 * np.pi, size)
y = np.sin(x) + np.random.random(size) * 0.2
data = {"x": x, "y": y}

Los paso a una función para medir el tiempo de ejecución y trazar el ajuste resultante:

def test_func(f, label):  # f: function handle to one of the smoothing methods
    start = time()
    for i in range(5):
        arr = f(data["y"], 20)
    print(f"{label:26s} - time: {time() - start:8.5f} ")
    plt.plot(data["x"], arr, "-", label=label)

Probé muchas funciones de suavizado diferentes. arr es la matriz de valores de y que se suavizará y span el parámetro de suavizado. Cuanto más bajo, mejor se acercará el ajuste a los datos originales, cuanto más alto, más suave será la curva resultante.

def smooth_data_convolve_my_average(arr, span):
    re = np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")

    # The "my_average" part: shrinks the averaging window on the side that 
    # reaches beyond the data, keeps the other side the same size as given 
    # by "span"
    re[0] = np.average(arr[:span])
    for i in range(1, span + 1):
        re[i] = np.average(arr[:i + span])
        re[-i] = np.average(arr[-i - span:])
    return re

def smooth_data_np_average(arr, span):  # my original, naive approach
    return [np.average(arr[val - span:val + span + 1]) for val in range(len(arr))]

def smooth_data_np_convolve(arr, span):
    return np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")

def smooth_data_np_cumsum_my_average(arr, span):
    cumsum_vec = np.cumsum(arr)
    moving_average = (cumsum_vec[2 * span:] - cumsum_vec[:-2 * span]) / (2 * span)

    # The "my_average" part again. Slightly different to before, because the
    # moving average from cumsum is shorter than the input and needs to be padded
    front, back = [np.average(arr[:span])], []
    for i in range(1, span):
        front.append(np.average(arr[:i + span]))
        back.insert(0, np.average(arr[-i - span:]))
    back.insert(0, np.average(arr[-2 * span:]))
    return np.concatenate((front, moving_average, back))

def smooth_data_lowess(arr, span):
    x = np.linspace(0, 1, len(arr))
    return sm.nonparametric.lowess(arr, x, frac=(5*span / len(arr)), return_sorted=False)

def smooth_data_kernel_regression(arr, span):
    # "span" smoothing parameter is ignored. If you know how to 
    # incorporate that with kernel regression, please comment below.
    kr = KernelReg(arr, np.linspace(0, 1, len(arr)), 'c')
    return kr.fit()[0]

def smooth_data_savgol_0(arr, span):  
    return savgol_filter(arr, span * 2 + 1, 0)

def smooth_data_savgol_1(arr, span):  
    return savgol_filter(arr, span * 2 + 1, 1)

def smooth_data_savgol_2(arr, span):  
    return savgol_filter(arr, span * 2 + 1, 2)

def smooth_data_fft(arr, span):  # the scaling of "span" is open to suggestions
    w = fftpack.rfft(arr)
    spectrum = w ** 2
    cutoff_idx = spectrum < (spectrum.max() * (1 - np.exp(-span / 2000)))
    w[cutoff_idx] = 0
    return fftpack.irfft(w)

Resultados

Velocidad

Tiempo de ejecución de más de 1000 elementos, probado en una lista de python, así como en una matriz numpy para contener los valores.

method              | python list | numpy array
--------------------|-------------|------------
kernel regression   | 23.93405 s  | 22.75967 s
lowess              |  0.61351 s  |  0.61524 s
numpy average       |  0.02485 s  |  0.02326 s
savgol 2            |  0.00186 s  |  0.00196 s
savgol 1            |  0.00157 s  |  0.00161 s
savgol 0            |  0.00155 s  |  0.00151 s
numpy convolve + me |  0.00121 s  |  0.00115 s
numpy cumsum + me   |  0.00114 s  |  0.00105 s
fft                 |  0.00021 s  |  0.00021 s
numpy convolve      |  0.00017 s  |  0.00015 s

Especialmente kernel regression es muy lento para calcular más de 1k elementos, lowess también falla cuando el conjunto de datos se vuelve mucho más grande. numpy convolve y fft son especialmente rápidos. No investigué el comportamiento del tiempo de ejecución (O(n)) con el aumento o la disminución del tamaño de la muestra.

Comportamiento de borde

Separaré esta parte en dos, para que la imagen sea comprensible.

Métodos basados ​​en Numpy + savgol 0:

Estos métodos calculan un promedio de los datos, el gráfico no se suaviza. Todos (a excepción de numpy.cumsum) dan como resultado el mismo gráfico cuando la ventana que se usa para calcular el promedio no toca el borde de los datos. la discrepancia a numpy.cumsum es muy probable que se deba a un error de ‘apagado por uno’ en el tamaño de la ventana.

Hay diferentes comportamientos de borde cuando el método tiene que trabajar con menos datos:

• savgol 0: continúa con una constante hasta el borde de los datos (savgol 1 y savgol 2 terminan con una recta y una parábola respectivamente)
• numpy average: se detiene cuando la ventana llega al lado izquierdo de los datos y llena esos lugares en la matriz con Nanmismo comportamiento que my_average método en el lado derecho
• numpy convolve: sigue los datos con bastante precisión. Sospecho que el tamaño de la ventana se reduce simétricamente cuando un lado de la ventana alcanza el borde de los datos
• my_average/me: mi propio método que implementé, porque no estaba satisfecho con los otros. Simplemente reduce la parte de la ventana que se extiende más allá de los datos hasta el borde de los datos, pero mantiene la ventana del otro lado del tamaño original dado con span

Métodos complicados:

Todos estos métodos terminan con un buen ajuste a los datos. savgol 1 termina con una línea, savgol 2 con una parábola.

Comportamiento de la curva

Mostrar el comportamiento de los diferentes métodos en medio de los datos.

Lo diferente savgol y average los filtros producen una línea áspera, lowess, fft y kernel regression producir un ajuste suave. lowess parece tomar atajos cuando cambian los datos.

Motivación

Tengo una Raspberry Pi que registra datos por diversión y la visualización resultó ser un pequeño desafío. Todos los puntos de datos, excepto el uso de RAM y el tráfico de Ethernet, solo se registran en pasos discretos y/o son inherentemente ruidosos. Por ejemplo, el sensor de temperatura solo emite grados enteros, pero difiere hasta en dos grados entre mediciones consecutivas. No se puede obtener información útil de un diagrama de dispersión de este tipo. Por lo tanto, para visualizar los datos, necesitaba algún método que no fuera demasiado costoso desde el punto de vista computacional y que produjera una media móvil. También quería un buen comportamiento en los bordes de los datos, ya que esto afecta especialmente la información más reciente cuando se miran los datos en vivo. me acomodé en el numpy convolve método con my_average para mejorar el comportamiento de los bordes.

Otra opción es utilizar KernelReg en modelos estadisticos:

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

Una definición clara de suavizado de una señal 1D de Libro de cocina SciPy le muestra cómo funciona.

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
    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
    xn=x+randn(len
    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()