如何以正确的方式平滑曲线？

Question

让我们假设我们有一个可能大约给出的数据集

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

因此,我们有20%的数据集变化.我的第一个想法是使用scipy的单变量函数函数,但问题是这不会很好地考虑小噪声.如果你考虑频率,背景远小于信号,所以只有截止的样条可能是一个想法,但这将涉及来回傅里叶变换,这可能导致不良行为.另一种方式是移动平均线,但这也需要正确选择延迟.

任何提示/书籍或链接如何解决这个问题？

Answer 1

我更喜欢Savitzky-Golay滤镜.它使用最小二乘法将数据的小窗口回归到多项式上,然后使用多项式估计窗口中心的点.最后,窗口向前移动一个数据点并重复该过程.这一直持续到每个点相对于其邻居进行了最佳调整.即使是来自非周期性和非线性源的噪声样本,它也能很好地工作.

这是一本完整的食谱示例.请参阅下面的代码,了解它的易用性.注意:我省略了用于定义savitzky_golay()函数的代码,因为您可以从上面链接的cookbook示例中逐字复制/粘贴它.

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

更新:我注意到我所连接的食谱示例已被删除.幸运的是,正如@dodohjk所指出的那样,Savitzky-Golay过滤器已被整合到SciPy库中.要使用SciPy源修改上述代码,请键入:

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

Answer 2

基于移动平均框(通过卷积)来平滑我使用的数据的快速而肮脏的方法:

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)

Answer 3

如果您对周期性信号(如您的示例)的"平滑"版本感兴趣,那么FFT是正确的方法.进行傅立叶变换并减去低贡献频率:

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)

即使您的信号不是完全周期性的,这也可以很好地减去白噪声.有许多类型的过滤器可供使用(高通,低通等),适当的过滤器取决于您所寻找的内容.

Answer 4

为您的数据拟合移动平均线可以消除噪音,请参阅此答案以了解如何执行此操作.

如果你想使用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()

最后,如果你知道信号的功能形式,你可以为你的数据拟合曲线,这可能是最好的事情.

Answer 5

这个问题已经得到了彻底的回答，所以我认为对所提议方法的运行时分析会很有趣（无论如何，这对我来说）。我还将在嘈杂数据集的中心和边缘查看方法的行为。

TL; 博士

                    | runtime in s | runtime in s
method              | python list  | numpy array
--------------------|--------------|------------
kernel regression   | 23.93405     | 22.75967 
lowess              |  0.61351     |  0.61524 
naive average       |  0.02485     |  0.02326 
others*             |  0.00150     |  0.00150 
fft                 |  0.00021     |  0.00021 
numpy convolve      |  0.00017     |  0.00015 

*savgol with different fit functions and some numpy methods

核回归的扩展性很差，Lowess 的速度要快一些，但两者都能产生平滑的曲线。Savgol 是速度的中间地带，可以产生跳跃和平滑的输出，具体取决于多项式的等级。FFT 非常快，但仅适用于周期性数据。

使用 numpy 的移动平均方法更快，但显然会生成一个带有步骤的图形。

设置

我以正弦曲线的形状生成了 1000 个数据点：

size = 1000
x = np.linspace(0, 4 * np.pi, size)
y = np.sin(x) + np.random.random(size) * 0.2
data = {"x": x, "y": y}

我将这些传递给一个函数来测量运行时间并绘制结果拟合：

def test_func(f, label):  # f: function handle to one of the smoothing methods
    start = time()
    for i in range(5):
        arr = f(data["y"], 20)
    print(f"{label:26s} - time: {time() - start:8.5f} ")
    plt.plot(data["x"], arr, "-", label=label)

我测试了许多不同的平滑功能。arr是要平滑的 y 值数组和span平滑参数。越低，拟合越接近原始数据，越高，结果曲线越平滑。

def smooth_data_convolve_my_average(arr, span):
    re = np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")

    # The "my_average" part: shrinks the averaging window on the side that 
    # reaches beyond the data, keeps the other side the same size as given 
    # by "span"
    re[0] = np.average(arr[:span])
    for i in range(1, span + 1):
        re[i] = np.average(arr[:i + span])
        re[-i] = np.average(arr[-i - span:])
    return re

def smooth_data_np_average(arr, span):  # my original, naive approach
    return [np.average(arr[val - span:val + span + 1]) for val in range(len(arr))]

def smooth_data_np_convolve(arr, span):
    return np.convolve(arr, np.ones(span * 2 + 1) / (span * 2 + 1), mode="same")

def smooth_data_np_cumsum_my_average(arr, span):
    cumsum_vec = np.cumsum(arr)
    moving_average = (cumsum_vec[2 * span:] - cumsum_vec[:-2 * span]) / (2 * span)

    # The "my_average" part again. Slightly different to before, because the
    # moving average from cumsum is shorter than the input and needs to be padded
    front, back = [np.average(arr[:span])], []
    for i in range(1, span):
        front.append(np.average(arr[:i + span]))
        back.insert(0, np.average(arr[-i - span:]))
    back.insert(0, np.average(arr[-2 * span:]))
    return np.concatenate((front, moving_average, back))

def smooth_data_lowess(arr, span):
    x = np.linspace(0, 1, len(arr))
    return sm.nonparametric.lowess(arr, x, frac=(5*span / len(arr)), return_sorted=False)

def smooth_data_kernel_regression(arr, span):
    # "span" smoothing parameter is ignored. If you know how to 
    # incorporate that with kernel regression, please comment below.
    kr = KernelReg(arr, np.linspace(0, 1, len(arr)), 'c')
    return kr.fit()[0]

def smooth_data_savgol_0(arr, span):  
    return savgol_filter(arr, span * 2 + 1, 0)

def smooth_data_savgol_1(arr, span):  
    return savgol_filter(arr, span * 2 + 1, 1)

def smooth_data_savgol_2(arr, span):  
    return savgol_filter(arr, span * 2 + 1, 2)

def smooth_data_fft(arr, span):  # the scaling of "span" is open to suggestions
    w = fftpack.rfft(arr)
    spectrum = w ** 2
    cutoff_idx = spectrum < (spectrum.max() * (1 - np.exp(-span / 2000)))
    w[cutoff_idx] = 0
    return fftpack.irfft(w)

结果

速度

运行时超过 1000 个元素，在 python 列表和 numpy 数组上测试以保存值。

method              | python list | numpy array
--------------------|-------------|------------
kernel regression   | 23.93405 s  | 22.75967 s
lowess              |  0.61351 s  |  0.61524 s
numpy average       |  0.02485 s  |  0.02326 s
savgol 2            |  0.00186 s  |  0.00196 s
savgol 1            |  0.00157 s  |  0.00161 s
savgol 0            |  0.00155 s  |  0.00151 s
numpy convolve + me |  0.00121 s  |  0.00115 s
numpy cumsum + me   |  0.00114 s  |  0.00105 s
fft                 |  0.00021 s  |  0.00021 s
numpy convolve      |  0.00017 s  |  0.00015 s

特别kernel regression是计算超过 1k 个元素lowess时非常慢，当数据集变得更大时也会失败。numpy convolve而且fft特别快。我没有研究增加或减少样本大小的运行时行为 (O(n))。

边缘行为

我将把这部分分成两部分，以保持图像易于理解。

基于 Numpy 的方法 + savgol 0：

这些方法计算数据的平均值，图形没有平滑。numpy.cumsum当用于计算平均值的窗口未触及数据边缘时，它们都（除了）生成相同的图形。与 的差异numpy.cumsum很可能是由于窗口大小中的“差一”错误造成的。

当该方法必须使用较少的数据时，会有不同的边缘行为：

• savgol 0：具有恒定的数据的边缘继续进行（savgol 1并savgol 2分别用线和抛物线结束）
• numpy average: 当窗口到达数据的左侧时停止并用 填充数组中的那些位置Nan，与my_average右侧的方法相同的行为
• numpy convolve: 非常准确地跟踪数据。我怀疑当窗口的一侧到达数据边缘时，窗口大小会对称减小
• my_average/ me: 我自己实现的方法，因为我对其他方法不满意。简单地将超出数据的窗口部分缩小到数据的边缘，但将窗口保持在另一侧的原始大小span

复杂的方法：

这些方法都以非常适合数据结束。savgol 1以一条直线结束，savgol 2一条抛物线。

曲线行为

在数据中间展示不同方法的行为。

不同的savgol和average过滤器产生一条粗线，lowess，fft并kernel regression产生平滑的配合。lowess当数据改变时，似乎偷工减料。

动机

我有一个用于娱乐的 Raspberry Pi 日志记录数据，事实证明可视化是一个小挑战。除 RAM 使用和以太网流量外，所有数据点仅以离散步骤和/或固有噪声记录。例如，温度传感器仅输出整度数，但在连续测量之间最多相差两度。从这样的散点图中无法获得有用的信息。因此，为了可视化数据，我需要一些计算成本不太高并产生移动平均线的方法。我还希望在数据边缘有良好的行为，因为这尤其会影响查看实时数据时的最新信息。我决定在numpy convolve与方法my_average，提高边缘的行为。

Answer 6

另一种选择是在statsmodel中使用KernelReg:

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

Answer 7

对于我的一个项目，我需要为时间序列建模创建间隔，为了使过程更加高效，我创建了tsmoothie：一个以矢量化方式进行时间序列平滑和异常值检测的 python 库。

它提供了不同的平滑算法以及计算间隔的可能性。

这里我使用的是一个，ConvolutionSmoother但你也可以测试其他的。

import numpy as np
import matplotlib.pyplot as plt
from tsmoothie.smoother import *

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

# operate smoothing
smoother = ConvolutionSmoother(window_len=5, window_type='ones')
smoother.smooth(y)

# generate intervals
low, up = smoother.get_intervals('sigma_interval', n_sigma=2)

# plot the smoothed timeseries with intervals
plt.figure(figsize=(11,6))
plt.plot(smoother.smooth_data[0], linewidth=3, color='blue')
plt.plot(smoother.data[0], '.k')
plt.fill_between(range(len(smoother.data[0])), low[0], up[0], alpha=0.3)

我还指出 tsmoothie 可以以矢量化的方式对多个时间序列进行平滑

Answer 8

看一下这个！对一维信号的平滑有一个明确的定义。

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