Как правильно сгладить кривую?

200

Предположим, у нас есть набор данных, который может быть дан примерно

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

Поэтому у нас есть вариация 20% набора данных. Моя первая идея состояла в том, чтобы использовать функцию Scivy UnivariateSpline, но проблема в том, что это не учитывает небольшой шум в хорошем смысле. Если вы рассматриваете частоты, фон намного меньше сигнала, поэтому сплайном только среза может быть идея, но это может включать в себя преобразование Фурье назад и вперед, что может привести к плохому поведению. Другим способом будет скользящее среднее, но для этого также потребуется правильный выбор задержки.

Любые советы / книги или ссылки, как решить эту проблему?

пример

varantir
источник
1
Ваш сигнал всегда будет синусоидальным, или вы использовали его только в качестве примера?
Марк Рэнсом
нет, у меня будут разные сигналы, даже в этом простом примере очевидно, что моих методов недостаточно
varantir
Фильтрация Калмана оптимальна для этого случая. И пакет pykalman python хорошего качества.
Тойня
Возможно, я расширю его до полного ответа, когда у меня будет немного больше времени, но один мощный метод регрессии, который еще не был упомянут, - это регрессия GP (Gaussian Process).
Ori5678

Ответы:

262

Я предпочитаю фильтр Савицкого-Голея . Он использует метод наименьших квадратов для регрессии небольшого окна ваших данных в полином, а затем использует полином для оценки точки в центре окна. Наконец, окно сдвигается вперед на одну точку данных, и процесс повторяется. Это продолжается до тех пор, пока каждая точка не будет оптимально отрегулирована относительно ее соседей. Он отлично работает даже с шумными выборками из непериодических и нелинейных источников.

Вот подробный пример поваренной книги . Посмотрите мой код ниже, чтобы понять, как легко им пользоваться. Примечание: я пропустил код для определения savitzky_golay()функции, потому что вы можете буквально скопировать / вставить ее из примера поваренной книги, который я привел выше.

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

оптимально сглаживать шумную синусоиду

ОБНОВЛЕНИЕ: до меня дошло, что пример поваренной книги, на который я ссылался, был удален. К счастью, фильтр Савицкого-Голея был включен в библиотеку SciPy , на что указывает @dodohjk . Чтобы адаптировать приведенный выше код с использованием исходного кода SciPy, введите:

from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
Дэвид Вюрц
источник
Я получил ошибку Traceback (последний вызов был последним): файл "hp.py", строка 79, в <module> ysm2 = savitzky_golay (y_data, 51,3) Файл "hp.py", строка 42, в savitzky_golay firstvals = y [0] - np.abs (y [1: half_window + 1] [:: - 1] - y [0])
март, Ho
14
Спасибо за введение фильтра Савицкого-Голея! Таким образом, в основном это похоже на обычный фильтр «Скользящее среднее», но вместо простого вычисления среднего, для каждой точки выполняется полиномиальный (обычно 2-й или 4-й порядок) подбор, и выбирается только «средняя» точка. Поскольку информация 2-го (или 4-го) порядка касается каждой точки, смещение, введенное в подходе «скользящее среднее» в локальных максимумах или минимумах, обходится. Действительно элегантный
np8
2
Просто хочу сказать спасибо за это, я схожу с ума, пытаясь выяснить вейвлет-разложения, чтобы получить сглаженные данные, и это намного приятнее.
Эльдар М.
5
Если х данные регулярно не разнесены вы можете применить фильтр к иксы , а также: savgol_filter((x, y), ...).
Тим Кейперс
127

Быстрый и грязный способ сглаживания данных, которые я использую, на основе блока скользящего среднего (путем свертки):

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)

введите описание изображения здесь

scrx2
источник
9
Это имеет несколько приятных преимуществ: (1) работает для любой функции, а не только для периодической, и (2) нет зависимостей или больших функций для копирования-вставки. Вы можете сделать это прямо сейчас с чистой Numpy. Кроме того, он не слишком грязный - это простейший случай некоторых из других методов, описанных выше (например, LOWESS, но ядро ​​имеет резкий интервал и, как Савицкий-Голей, но степень полинома равна нулю).
Джим Пиварски
2
единственная проблема со скользящей средней состоит в том, что она отстает от данных. Вы можете увидеть это наиболее очевидно в конце, где больше точек сверху и меньше снизу, но зеленая кривая в настоящее время ниже среднего, потому что оконная функция должна двигаться вперед, чтобы учесть это.
Нуреттин
И это не работает на массиве nd, только 1d. scipy.ndimage.filters.convolve1d()позволяет указать ось nd-массива для фильтрации. Но я думаю, что оба страдают от некоторых проблем в замаскированных ценностях.
Джейсон
1
@nurettin Я думаю, что вы описываете краевые эффекты. В общем, до тех пор, пока ядро ​​свертки способно покрыть свою протяженность в сигнале, оно, как вы говорите, не «отстает». В конце, однако, нет значений, превышающих 6, для включения в среднее значение, поэтому используется только «левая» часть ядра. Краевые эффекты присутствуют в каждом ядре сглаживания и должны обрабатываться отдельно.
Джон
4
@nurettin Нет, я пытался разъяснить другим читателям, что ваш комментарий «единственная проблема с скользящей средней - это то, что она отстает от данных», вводит в заблуждение. Любой метод оконного фильтра страдает этой проблемой, а не только скользящим средним. Савицкий-Голай тоже страдает этой проблемой. Таким образом, ваше утверждение «То, что я описываю, это то, что savitzky_golay решает с помощью оценки», просто неверно. Любой метод сглаживания требует способа обработки краев, который не зависит от самого метода сглаживания.
Джон
79

Если вы заинтересованы в «гладкой» версии сигнала, которая является периодической (как ваш пример), тогда БПФ - правильный путь. Возьмите преобразование Фурье и вычтите низкочастотные частоты:

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)

введите описание изображения здесь

Даже если ваш сигнал не является полностью периодическим, это позволит отлично вычесть белый шум. Существует много типов фильтров (высокочастотный, низкочастотный и т. Д.), Подходящий из которых зависит от того, что вы ищете.

Увлеченные
источник
Какой участок для какой переменной? Я пытаюсь сгладить координаты для теннисного мяча в ралли, т.е.
убери все скачки,
44

Подгонка скользящего среднего к вашим данным сгладит шум, посмотрите этот ответ, чтобы узнать, как это сделать.

Если вы хотите использовать LOWESS, чтобы соответствовать вашим данным (это похоже на скользящее среднее, но более изощренно), вы можете сделать это с помощью библиотеки 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()

Наконец, если вы знаете функциональную форму вашего сигнала, вы можете подогнать кривую к вашим данным, что, вероятно, будет лучшим решением.

markmuetz
источник
Если бы только это было loessреализовано.
scrutari
18

Другой вариант - использовать KernelReg в 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()
Зичен Ван
источник
7

Проверь это! Существует четкое определение сглаживания 1D сигнала.

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

Клавиши быстрого доступа:

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
источник
3
Ссылка на решение приветствуется, но, пожалуйста, убедитесь, что ваш ответ будет полезен без него: добавьте контекст вокруг ссылки, чтобы ваши коллеги-пользователи поняли, что это такое и почему, а затем процитируйте наиболее релевантную часть страницы, которую вы Повторная ссылка на случай, если целевая страница недоступна. Ответы, которые немного больше, чем ссылка, могут быть удалены.
Шри
-4

Если вы строите график временных рядов и используете mtplotlib для рисования графиков, используйте метод медианы для сглаживания графика.

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

где timeseriesваш набор данных передается вы можете изменить windowsizeдля более сглаживания.

Паван Пурохит
источник