logo标签列表

如何为数据集平滑曲线

pythonnumpyscipysignal-processingsmoothing
294
假设我们有一个数据集，可能大致如下所示
import numpy as np
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2

因此，我们的数据集有20%的变化。我的第一个想法是使用scipy的UnivariateSpline函数，但问题是它不能很好地考虑到小噪声。如果考虑频率，背景远小于信号，所以只对截断部分进行样条插值可能是个主意，但这将涉及前后傅里叶变换，可能导致不良行为。 另一种方法是移动平均，但这也需要正确选择延迟。
有没有什么提示/书籍或链接可以解决这个问题？

example

12个回答
389

我更喜欢使用 Savitzky-Golay 滤波器。它利用最小二乘法将数据的一个小窗口回归到一个多项式，然后使用该多项式来估计窗口中心的点。最后，窗口向前移动一个数据点，重复该过程。这一过程会持续进行，直到每个点都相对于其邻居进行了最优调整。即使是从非周期性和非线性源获得的嘈杂样本，它也能很好地工作。

这里有一个详尽的示例。查看下面的代码以了解使用起来有多简单。注意：我省略了定义 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()

平滑处理噪声正弦波的最佳方法

更新：我注意到我链接的示例已被删除。幸运的是，Savitzky-Golay 滤波器已经被整合进入SciPy库，如@dodohjk指出（感谢@bicarlsen提供的更新链接）。 要使用SciPy源代码来调整上面的代码，请输入：

from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
2
我得到了错误回溯（最近的调用最先）： 在“hp.py”中第79行，ysm2 = savitzky_golay（y_data，51，3） 在“hp.py”中第42行，firstvals = y [0] - np.abs（y [1：half_window + 1] [::-1] - y [0]） - March Ho
11
如果 x 数据不是等间距的，您可能希望将该滤波器应用于 x 值：savgol_filter((x, y), ...)。 - Tim Kuipers
2
说一个程序可以与“非线性源”一起工作是什么意思？什么是“非线性源”？ - DanielSank
2
@TimKuipers 我尝试了这个，但是出现了错误，因为现在x参数只有大小2（scipy函数似乎没有“深入”查看这实际上是一个大小为m的数组元组，对于m个数据点中的每个数组） - Chris K
2
scipy.signal#savgol_filter的链接已经失效，但我相信这是正确的链接：https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html - bicarlsen
162

编辑：查看此答案。使用np.cumsum比使用np.convolve更快。

我用一种基于移动平均箱（通过卷积）的快速而简单的方法来平滑数据：

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)

输入图像描述

1
这在nd数组上不起作用，只能在1d上使用。scipy.ndimage.filters.convolve1d()允许您指定一个nd数组的轴来进行过滤。但我认为两者都存在一些掩码值问题。 - Jason
6
我在数组的开头和结尾处出现奇怪的边缘效应（第一个和最后一个值约为其他值的一半）。 - Chris
92

如果你对一个周期性信号（像你的例子）有兴趣获得“平滑”版本的话，那么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)

enter image description here

即使您的信号不完全周期性，这也会很好地消除白噪声。有许多类型的滤波器可供使用（高通、低通等），适当的滤波器取决于您要查找的内容。

1
哪个图是哪个变量？我正在尝试平滑网球比赛中球的坐标，即从我的图表中删除所有看起来像小抛物线的弹跳。 - mLstudent33
59

这个问题已经得到了彻底的回答，因此我认为对所提出方法的运行时分析会很有兴趣（对我来说确实是这样）。 我还将研究这些方法在嘈杂数据集的中心和边缘的行为。

简而言之

                    | 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

Kernel回归的缩放效果不佳，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)

结果

速度

测试了一个python列表和一个numpy数组来保存值，对超过1000个元素的运行时间进行了测试。

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

特别是核回归在计算超过1k个元素时非常缓慢，当数据集变得更大时，lowess也会失败。尤其是numpy convolve和fft非常快速。我没有研究随着样本大小增加或减少的运行时行为（O(n)）。

边缘行为

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

Numpy基于方法+ savgol 0：

numpy based methods的边缘行为

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

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

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

复杂的方法： 复杂方法的边缘行为

这些方法都以数据的良好拟合结束。savgol 1以直线结束，savgol 2以抛物线结束。

曲线行为

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

不同方法的曲线行为

不同的savgol和average滤波器产生粗糙的直线，lowess、fft和kernel regression产生平滑的拟合。当数据变化时，lowess似乎切角。

动机

我有一个树莓派记录数据以供娱乐，可视化证明是一个小挑战。除了RAM使用情况和以太网流量外，所有数据

这是一个非常详细的回复 - 谢谢！我想通过可视化函数来了解Convolve平滑与my_average的作用...将尝试在matplotlib上构建它... - Pushkar Sheth
58

对你的数据进行移动平均拟合可以平滑噪声，可查看此答案了解如何操作。

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

最后，如果您知道信号的功能形式，您可以将曲线拟合到数据上，这可能是最好的选择。

如果只有loess被实现了。 - scrutari
23

另一种选择是使用statsmodels中的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()
KernalReg不能平滑曲线。 - ERIC
KernalReg确实平滑了曲线。它是我问题的首选！ - Syrtis Major
9

从SciPy Cookbook中清晰地定义了一维信号平滑的概念，它展示了其工作原理。

快捷方式：

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()
3
欢迎提供解决方案的链接，但请确保您的回答在没有链接的情况下也有用：添加链接周围的上下文，以便其他用户了解它是什么以及为什么存在，然后引用您链接页面中与问题最相关的部分，以防目标页面无法访问。只有链接而没有更多内容的答案可能会被删除。 - 4b0
6
为了我的项目，我需要为时间序列建立间隔，为了使该过程更加高效，我创建了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)

enter image description here

我还要指出，tsmoothie可以以向量化的方式对多个时间序列进行平滑处理。
2

scipy.ndimage 中有一个简单的函数，对我来说也很有效：

from scipy.ndimage import uniform_filter1d

y_smooth = uniform_filter1d(y,size=15)

enter image description here

0
你也可以使用这个：
def smooth(scalars, weight = 0.8):  # Weight between 0 and 1
        return [scalars[i] * weight + (1 - weight) * scalars[i+1] for i in range(len(scalars)) if i < len(scalars)-1]
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接