Jak wygładzić krzywą we właściwy sposób?

Załóżmy, że mamy zestaw danych, który może być podany w przybliżeniu przez

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

Dlatego mamy zmianę o 20% zbioru danych. Moim pierwszym pomysłem było użycie funkcji scipy UnivariateSpline, ale problem polega na tym, że nie uwzględnia to właściwie małego szumu. Jeśli weźmie się pod uwagę częstotliwości, tło jest znacznie mniejsze niż sygnał, więc pomysłem może być tylko splajn odcięcia, ale wymagałoby to transformacji Fouriera w przód iw tył, co może skutkować złym zachowaniem. Innym sposobem byłaby średnia ruchoma, ale wymagałoby to również właściwego wyboru opóźnienia.

Wszelkie wskazówki / książki lub linki, jak rozwiązać ten problem?

Wolę filtr Savitzky-Golay . Używa najmniejszych kwadratów, aby regresować małe okno danych do wielomianu, a następnie używa wielomianu do oszacowania punktu na środku okna. Na koniec okno zostaje przesunięte o jeden punkt danych do przodu i proces się powtarza. Trwa to do momentu, aż każdy punkt zostanie optymalnie dostosowany w stosunku do swoich sąsiadów. Działa świetnie nawet w przypadku głośnych próbek ze źródeł nieokresowych i nieliniowych.

Oto dokładny przykład książki kucharskiej . Zobacz mój kod poniżej, aby dowiedzieć się, jak łatwo jest go używać. Uwaga: Pominąłem kod definiujący savitzky_golay()funkcję, ponieważ możesz dosłownie skopiować / wkleić go z przykładu książki kucharskiej, który podałem powyżej.

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

AKTUALIZACJA: Zwróciłem uwagę, że przykład książki kucharskiej, z którym się łączyłem, został usunięty. Na szczęście filtr Savitzky-Golay został włączony do biblioteki SciPy , jak wskazał @dodohjk . Aby dostosować powyższy kod za pomocą źródła SciPy, wpisz:

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

Szybki i brudny sposób na wygładzanie danych, którego używam, w oparciu o średnią ruchomą (przez splot):

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)

Jeśli interesuje Cię „gładka” wersja sygnału, który jest okresowy (jak twój przykład), to FFT jest właściwą drogą. Weź transformatę Fouriera i odejmij niskie częstotliwości:

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)

Nawet jeśli twój sygnał nie jest całkowicie okresowy, to doskonale odciąży biały szum. Istnieje wiele rodzajów filtrów do użycia (górnoprzepustowy, dolnoprzepustowy itp.), Odpowiedni jest zależny od tego, czego szukasz.

Dopasowanie średniej ruchomej do danych zmniejszy hałas, zapoznaj się z tą odpowiedzią, jak to zrobić.

Jeśli chcesz użyć opcji LOWESS, aby dopasować swoje dane (jest podobny do średniej ruchomej, ale bardziej wyrafinowany), możesz to zrobić za pomocą biblioteki 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()

Wreszcie, jeśli znasz funkcjonalną formę sygnału, możesz dopasować krzywą do swoich danych, co prawdopodobnie byłoby najlepszą rzeczą do zrobienia.

Inną opcją jest użycie KernelReg w 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()

Sprawdź to! Istnieje jasna definicja wygładzania sygnału 1D.

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

Skrót:

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

Jeśli drukujesz wykres szeregów czasowych i użyłeś mtplotlib do rysowania wykresów, użyj metody mediany, aby wygładzić wykres

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

gdzie timeseriesjest przekazywany zestaw danych, który możesz zmienić, windowsizeaby uzyskać bardziej wygładzenie.