wake-up-neo.net

Wie glättet man eine Kurve richtig?

Nehmen wir an, wir haben einen Datensatz, der ungefähr durch angegeben werden kann

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

Daher haben wir eine Abweichung von 20% des Datensatzes. Meine erste Idee war, die UnivariateSpline-Funktion von scipy zu verwenden, aber das Problem ist, dass dies das kleine Rauschen nicht gut berücksichtigt. Wenn Sie die Frequenzen in Betracht ziehen, ist der Hintergrund viel kleiner als das Signal. Daher kann ein Spline nur der Cutoff-Funktion eine Idee sein. Dies würde jedoch eine hin- und hergehende Fourier-Transformation bedeuten, die zu einem schlechten Verhalten führen könnte wäre ein gleitender Durchschnitt, aber dies würde auch die richtige Wahl der Verzögerung erfordern.

Irgendwelche Tipps/Bücher oder Links, um dieses Problem anzugehen?

example

131
varantir

Ich bevorzuge einen Savitzky-Golay-Filter . Es verwendet kleinste Quadrate, um ein kleines Fenster Ihrer Daten auf ein Polynom zurückzuführen, und dann das Polynom, um den Punkt in der Mitte des Fensters zu schätzen. Schließlich wird das Fenster um einen Datenpunkt nach vorne verschoben und der Vorgang wiederholt sich. Dies setzt sich fort, bis jeder Punkt im Verhältnis zu seinen Nachbarn optimal eingestellt ist. Es funktioniert auch mit geräuschvollen Samples aus nicht periodischen und nichtlinearen Quellen.

Hier ist ein ausführliches Kochbuch-Beispiel . Sehen Sie sich meinen Code unten an, um eine Vorstellung davon zu bekommen, wie einfach es zu benutzen ist. Hinweis: Ich habe den Code zur Definition der Funktion savitzky_golay() weggelassen, da Sie ihn aus dem oben verlinkten Kochbuchbeispiel kopieren und einfügen können.

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

optimally smoothing a noisy sinusoid

UPDATE: Es ist mir aufgefallen, dass das Kochbuchbeispiel, mit dem ich verlinkt habe, entfernt wurde. Glücklicherweise wurde der Savitzky-Golay-Filter in in die SciPy-Bibliothek eingebaut, wie unter @dodohjk ..__ angegeben. Um den obigen Code mithilfe der SciPy-Quelle anzupassen, geben Sie Folgendes ein:

from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
170
David Wurtz

Eine schnelle und schmutzige Methode zum Glätten von Daten, die ich verwende, basierend auf einem gleitenden Durchschnittswert (durch Faltung):

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)

 enter image description here

87
scrx2

Wenn Sie an einer "glatten" Version eines Signals interessiert sind (wie in Ihrem Beispiel), dann ist eine FFT der richtige Weg. Nehmen Sie die Fourier-Transformation und subtrahieren Sie die Frequenzen mit niedrigem Beitrag:

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)

enter image description here

Selbst wenn Ihr Signal nicht vollständig periodisch ist, können Sie damit weißes Rauschen herausziehen. Es gibt viele Arten von Filtern (Hochpass, Tiefpass usw.), der passende hängt von dem ab, was Sie suchen.

69
Hooked

Wenn Sie einen gleitenden Durchschnitt an Ihre Daten anpassen, wird das Rauschen ausgeglichen. Siehe diese Antwort , um zu erfahren, wie Sie dies tun können.

Wenn Sie LOWESS verwenden möchten, um Ihre Daten anzupassen (dies ähnelt einem gleitenden Durchschnitt, ist jedoch komplexer), können Sie dies mit der Bibliothek statsmodels tun:

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

Wenn Sie die Funktionsform Ihres Signals kennen, können Sie schließlich eine Kurve an Ihre Daten anpassen, was wahrscheinlich das Beste ist.

37
markmuetz

Eine andere Option ist die Verwendung von KernelReg in 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()
13
Zichen Wang

Schau dir das an! Die Glättung eines 1D-Signals ist klar definiert.

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

Abkürzung:

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