python - help - scipy reference

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

Question

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 utilizar la función UnivariateSpline de scipy, pero el problema es que esto no considera el ruido pequeño en una buena forma. Si se consideran las frecuencias, el fondo es mucho más pequeño que la señal, por lo que una spline solo del punto de corte podría ser una idea, pero eso implicaría una transformación fourier de ida y vuelta, lo que podría resultar en un mal comportamiento. Otra forma sería una media móvil, pero esto también necesitaría la elección correcta de la demora.

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

Answer 1

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

Si desea utilizar LOWESS para ajustarse a sus datos (es similar a una media móvil pero más sofisticada), puede hacerlo utilizando la biblioteca 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 sea lo mejor que puede hacer.

Answer 2

Otra opción es usar KernelReg en statsmodel :

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

Answer 3

Prefiero un filtro Savitzky-Golay . Utiliza mínimos cuadrados para hacer retroceder una pequeña ventana de tus 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 en un punto de datos y el proceso se repite. Esto continúa hasta que cada punto se haya ajustado de forma ó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 buen ejemplo de libro de cocina . Vea mi código a continuación para tener una idea de lo fácil que es usar. Nota: savitzky_golay() el código para definir la función savitzky_golay() porque literalmente puede copiarlo / pegarlo del ejemplo del libro de cocina que he vinculado 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()

ACTUALIZACIÓN: Me ha llamado la atención que el ejemplo del libro de cocina al que me he vinculado ha sido eliminado. Afortunadamente, parece que el filtro Savitzky-Golay se ha incorporado a la biblioteca SciPy , como lo señala @dodohjk .

Answer 4

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 a seguir. 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 restar ruido blanco. Hay muchos tipos de filtros para usar (paso alto, paso bajo, etc.), el adecuado depende de lo que esté buscando.

Answer 5

Una manera rápida y sucia de suavizar los datos que uso, en función de 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)