Como suavizar uma curva da maneira certa?

200

Vamos supor que temos um conjunto de dados que pode ser fornecido aproximadamente por

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

Portanto, temos uma variação de 20% do conjunto de dados. Minha primeira ideia foi usar a função UnivariateSpline do scipy, mas o problema é que isso não considera o pequeno ruído de uma maneira boa. Se você considerar as frequências, o fundo é muito menor que o sinal, portanto, uma spline apenas do ponto de corte pode ser uma ideia, mas isso envolveria uma transformação de quatro para a frente e para trás, o que poderia resultar em mau comportamento. Outra maneira seria uma média móvel, mas isso também precisaria da escolha certa do atraso.

Alguma dica / livros ou links sobre como lidar com esse problema?

exemplo

varantir
fonte
1
Seu sinal sempre será uma onda senoidal ou você estava usando isso apenas como exemplo?
precisa
não, eu terei sinais diferentes, mesmo neste exemplo fácil, é óbvio que meus métodos não são suficientes
varantir
A filtragem kalman é ideal para este caso. E o pacote pykalman python é de boa qualidade.
toine 23/02
Talvez eu o expanda para uma resposta completa quando tiver mais tempo, mas o método de regressão poderoso que ainda não foi mencionado é a regressão GP (Processo Gaussiano).
Ori5678 24/01

Respostas:

262

Eu prefiro um filtro Savitzky-Golay . Ele usa mínimos quadrados para regredir uma pequena janela dos seus dados em um polinômio e, em seguida, usa o polinômio para estimar o ponto no centro da janela. Finalmente, a janela é deslocada para frente em um ponto de dados e o processo se repete. Isso continua até que cada ponto seja ajustado de maneira ideal em relação aos seus vizinhos. Funciona muito bem mesmo com amostras barulhentas de fontes não periódicas e não lineares.

Aqui está um exemplo completo do livro de receitas . Veja meu código abaixo para ter uma idéia de como é fácil usar. Nota: Eu deixei de fora o código para definir a savitzky_golay()função porque você pode literalmente copiá-lo / colá-lo do exemplo do livro de receitas que eu vinculei acima.

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

otimizar suavemente um sinusóide barulhento

ATUALIZAÇÃO: Chegou ao meu conhecimento que o exemplo do livro de receitas ao qual vinculei foi retirado. Felizmente, o filtro Savitzky-Golay foi incorporado à biblioteca SciPy , como apontado por @dodohjk . Para adaptar o código acima usando a fonte SciPy, digite:

from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
David Wurtz
fonte
Eu recebi o erro Traceback (última chamada mais recente): arquivo "hp.py", linha 79, em <module> ysm2 = savitzky_golay (y_data, 51,3) Arquivo "hp.py", linha 42, em savitzky_golay firstvals = y [0] - np.abs (y [1: half_window + 1] [:: - 1] - y [0])
março Ho
14
Obrigado por apresentar o filtro Savitzky-Golay! Então, basicamente, é como um filtro "Média móvel" regular, mas em vez de apenas calcular a média, um ajuste polinomial (geralmente de 2ª ou 4ª ordem) é feito para cada ponto, e apenas o ponto "médio" é escolhido. Como as informações de 2º (ou 4º) pedidos estão preocupadas em todos os pontos, o viés introduzido na abordagem de "média móvel" em máximos ou mínimos locais é contornado. Realmente elegante.
Np8 26/03
2
Só quero dizer obrigado por isso, eu fiquei louca tentando descobrir decomposições de wavelets para obter dados suavizados, e isso é muito melhor.
Eldar M.
5
Se os dados x não é espaçadas regularmente você pode querer aplicar o filtro das x, bem como: savgol_filter((x, y), ...).
precisa saber é o seguinte
127

Uma maneira rápida e suja de suavizar os dados que eu uso, com base em uma caixa de média móvel (por convolução):

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)

insira a descrição da imagem aqui

scrx2
fonte
9
Isso tem algumas vantagens interessantes: (1) funciona para qualquer função, não apenas periódica, e (2) nenhuma dependência ou função grande para copiar e colar. Você pode fazer isso imediatamente com o Numpy puro. Além disso, não é muito sujo - é um caso mais simples de alguns dos outros métodos descritos acima (como LOWESS, mas o kernel é um intervalo nítido e como Savitzky-Golay, mas o grau polinomial é zero).
22615 Jim Jim Pivarski
2
o único problema com a média móvel é que fica atrás dos dados. Você pode ver isso mais aparentemente no final, onde há mais pontos na parte superior e menos na parte inferior, mas atualmente a curva verde está abaixo da média porque a função da janela precisa avançar para levar em consideração.
Nurettin
E isso não funciona na matriz nd, apenas 1d. scipy.ndimage.filters.convolve1d()permite especificar um eixo de um nd-array para fazer a filtragem. Mas acho que ambos sofrem de alguns problemas nos valores mascarados.
Jason
1
@ Nurettin Acho que o que você está descrevendo são efeitos de borda. Em geral, desde que o kernel de convolução seja capaz de cobrir sua extensão dentro do sinal, ele não ficará "atrasado", como você diz. No final, no entanto, não há valores além de 6 para incluir na média, portanto, apenas a parte "esquerda" do kernel está sendo usada. Os efeitos de borda estão presentes em todos os kernel de suavização e devem ser tratados separadamente.
Jon
4
@ Nurettin Não, eu estava tentando esclarecer para outras pessoas que estão lendo isso que seu comentário "o único problema com a média móvel é que fica atrás dos dados" é enganoso. Qualquer método de filtro de janela sofre esse problema, não apenas a média móvel. Savitzky-golay também sofre desse problema. Portanto, sua afirmação "O que estou descrevendo é o que savitzky_golay resolve por estimativa" está errada. Qualquer método de suavização requer uma maneira de lidar com arestas independente do próprio método de suavização.
Jon
79

Se você estiver interessado em uma versão "suave" de um sinal periódico (como o seu exemplo), uma FFT é o caminho certo a seguir. Faça a transformação de Fourier e subtraia as frequências de baixa contribuição:

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)

insira a descrição da imagem aqui

Mesmo que seu sinal não seja completamente periódico, isso fará um ótimo trabalho em subtrair o ruído branco. Existem muitos tipos de filtros a serem usados ​​(passa-alto, passa-baixo, etc ...), o apropriado depende do que você está procurando.

Hooked
fonte
Qual plot é para qual variável? Estou tentando suavizar as coordenadas para a bola de tênis em um comício, ou seja. tire todas as rejeições que pareçam pequenas parábolas no meu enredo
mLstudent33 03/03
44

Ajustar uma média móvel aos seus dados reduziria o ruído; veja esta resposta para saber como fazer isso.

Se você deseja usar o LOWESS para ajustar seus dados (é semelhante a uma média móvel, mas mais sofisticada), você pode fazer isso usando a 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()

Por fim, se você conhece a forma funcional do seu sinal, pode ajustar uma curva aos seus dados, o que provavelmente seria a melhor coisa a fazer.

markmuetz
fonte
Se apenas o tivesse loessimplementado.
Scrutari
18

Outra opção é usar o KernelReg no 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
fonte
7

Veja isso! Há uma definição clara de suavização de um sinal 1D.

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

Atalho:

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
fonte
3
Um link para uma solução é bem-vindo, mas garanta que sua resposta seja útil sem ela: adicione contexto ao link para que seus colegas usuários tenham uma idéia do que é e por que está lá, depois cite a parte mais relevante da página que você ' reencaminhando para o caso de a página de destino não estar disponível. Respostas que são pouco mais que um link podem ser excluídas.
Shree
-4

Se você estiver plotando um gráfico de séries temporais e se tiver usado o mtplotlib para desenhar gráficos, use o método mediano para suavizar o gráfico

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

onde timeseriesseu conjunto de dados é passado, você pode alterar windowsizepara obter mais suavização.

Pavan Purohit
fonte