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

200

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?

przykład

varantir
źródło
1
Czy twój sygnał zawsze będzie falą sinusoidalną, czy używałeś go tylko dla przykładu?
Mark Ransom,
nie, będę miał różne sygnały, nawet w tym prostym przykładzie jest oczywiste, że moje metody nie są wystarczające
varantir
Filtrowanie Kalmana jest optymalne w tym przypadku. Pakiet pykalman python jest dobrej jakości.
toine
Może rozwinę ją do pełnej odpowiedzi, gdy będę miał trochę więcej czasu, ale jedyną potężną metodą regresji, o której jeszcze nie wspomniano, jest regresja GP (Proces Gaussa).
Ori5678

Odpowiedzi:

261

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

optymalnie wygładza głośny sinusoidę

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
David Wurtz
źródło
Mam błąd Traceback (ostatnie ostatnie połączenie): Plik „hp.py”, wiersz 79, w <module> ysm2 = savitzky_golay (y_data, 51,3) Plik „hp.py”, wiersz 42, w savitzky_golay firstvals = y [0] - np.abs (y [1: half_window + 1] [:: - 1] - y [0])
March Ho
14
Dziękujemy za wprowadzenie filtra Savitzky-Golay! Zasadniczo jest to więc zwykły filtr „średniej ruchomej”, ale zamiast tylko obliczać średnią, dopasowuje się wielomian (zwykle 2. lub 4. rzędu) dla każdego punktu i wybierany jest tylko punkt „środkowy”. Ponieważ informacje o drugim (lub czwartym) rzędzie dotyczą każdego punktu, odchylenie wprowadzone w podejściu „średniej ruchomej” przy lokalnych maksimach lub minimach jest obchodzone. Naprawdę elegancki.
np8
2
Chcę tylko podziękować za to, że zwariowałem, próbując znaleźć rozkład falek, aby uzyskać wygładzone dane, a to jest o wiele ładniejsze.
Eldar M.
5
Jeśli dane x nie jest regularnie rozmieszczone może chcesz zastosować filtr do x, jak również: savgol_filter((x, y), ...).
Tim Kuipers,
127

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)

wprowadź opis zdjęcia tutaj

scrx2
źródło
9
Ma to kilka zalet: (1) działa dla dowolnej funkcji, nie tylko okresowej, i (2) nie ma zależności ani dużych funkcji do kopiowania i wklejania. Możesz to zrobić od razu za pomocą czystej Numpy. Ponadto, nie jest zbyt brudny - to najprostszy przypadek niektórych innych metod opisanych powyżej (takich jak LOWESS, ale jądro ma ostry interwał i jak Savitzky-Golay, ale stopień wielomianu wynosi zero).
Jim Pivarski
2
jedynym problemem związanym ze średnią ruchomą jest to, że pozostaje ona w tyle za danymi. Widać to najwyraźniej na końcu, gdzie na górze jest więcej kropek, a na dole mniej, ale zielona krzywa jest obecnie poniżej średniej, ponieważ funkcja okna musi przejść do przodu, aby je uwzględnić.
nurettin
I to nie działa na nd array, tylko 1d. scipy.ndimage.filters.convolve1d()pozwala określić oś nd-tablicy, aby wykonać filtrowanie. Ale myślę, że oboje cierpią na pewne problemy związane z zamaskowanymi wartościami.
Jason
1
@ nurettin Myślę, że opisujesz efekty krawędzi. Ogólnie, dopóki jądro splotu jest w stanie pokryć swój zasięg w sygnale, nie „opóźnia się”, jak mówisz. Na koniec jednak nie ma wartości wyższych niż 6, które należy uwzględnić w średniej, więc używana jest tylko „lewa” część jądra. Efekty krawędzi są obecne w każdym jądrze wygładzającym i muszą być obsługiwane osobno.
Jon
4
@nurettin Nie, próbowałem wyjaśnić innym czytającym to, że twój komentarz „jedynym problemem z ruchomą średnią jest to, że pozostaje w tyle za danymi” jest mylący. Każda metoda filtrowania okien ma ten problem, nie tylko ruchomą średnią. Savitzky-golay również cierpi na ten problem. Tak więc twoje stwierdzenie „To, co opisuję, jest tym, co savitzky_golay rozwiązuje przez oszacowanie” jest po prostu błędne. Każda z metod wygładzania wymaga sposobu obsługi krawędzi, który jest niezależny od samej metody wygładzania.
Jon
79

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)

wprowadź opis zdjęcia tutaj

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.

Haczykowaty
źródło
Który wykres jest dla jakiej zmiennej? Próbuję wygładzić współrzędne piłki tenisowej w rajdzie, tj.
usuń wszystkie odbicia,
44

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.

markmuetz
źródło
Gdyby tylko loesswdrożył.
scrutari
18

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()
Zichen Wang
źródło
7

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()
IPhysResearch
źródło
3
Link do rozwiązania jest mile widziany, ale upewnij się, że Twoja odpowiedź jest przydatna bez niego: dodaj kontekst wokół linku, aby inni użytkownicy mieli pojęcie, co to jest i dlaczego tam jest, a następnie zacytuj najbardziej odpowiednią część strony, którą „ ponowne linkowanie na wypadek, gdyby strona docelowa była niedostępna. Odpowiedzi, które są niewiele więcej niż link, mogą zostać usunięte.
Shree
-4

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.

Pavan Purohit
źródło