小开

如果你对周期信号的“平滑”版本感兴趣(就像你的例子)，那么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

即使你的信号不是完全周期性的，这也能很好地去除白噪声。有许多类型的过滤器可以使用(高通，低通，等等…)，合适的一个取决于你正在寻找什么。

小开

为你的数据拟合一个移动平均线可以消除噪音，请参阅这个答案了解如何做到这一点。

如果你想使用洛斯来拟合你的数据(它类似于移动平均，但更复杂)，你可以使用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()

最后，如果你知道信号的函数形式，你就可以为你的数据拟合曲线，这可能是最好的办法。

小开

我更喜欢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

小开

另一个选项是在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()

小开

编辑:看这答案。使用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)

小开

一个来自SciPy食谱的1D信号平滑的清晰定义向你展示了它是如何工作的。

快捷方式:

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

小开

如果你正在绘制时间序列图，如果你已经使用mtplotlib来绘制图形，那么使用中值方法平滑图

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

其中timeseries是你传递的数据集，你可以改变windowsize以进行更平滑。

小开

这个问题已经得到了彻底的回答，所以我认为对所提议的方法进行运行时分析会很有兴趣(至少对我来说是这样)。我还将查看噪声数据集中心和边缘的方法的行为。

博士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

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)

结果

速度

运行时超过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的差异很可能是由于窗口大小的“off by 1”错误。

当方法必须处理更少的数据时，有不同的边缘行为:

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

< p >复杂的方法:

< / p >

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

曲线的行为

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

不同的savgol和average过滤器产生粗线，lowess， fft和kernel regression产生平滑拟合。当数据发生变化时，lowess似乎可以走捷径。

动机

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

小开

对于我的一个项目，我需要为时间序列建模创建间隔，为了使过程更有效，我创建了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可以用向量化的方式对多个时间序列进行平滑

小开

使用移动平均线，一种快速的方法(也适用于非双射函数)是

def smoothen(x, winsize=5):
return np.array(pd.Series(x).rolling(winsize).mean())[winsize-1:]

此代码基于https://towardsdatascience.com/data-smoothing-for-data-science-visualization-the-goldilocks-trio-part-1-867765050615。文中还讨论了更先进的解决方案。

小开

你也可以用这个:

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]

小开

在scipy.ndimage中有一个简单的函数对我来说也很有用:

from scipy.ndimage import uniform_filter1d


y_smooth = uniform_filter1d(y,size=15)

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

博士TL;

设置

结果

速度

边缘的行为

曲线的行为

动机