平滑曲线的算法，严格保持在原曲线之上

Question

我有一个任意的输入曲线，以numpy数组的形式给出。 我想创建一个平滑版本，类似于滚动均值，但严格大于原始版本并且严格平滑。 我可以使用滚动均值，但是如果输入曲线具有负峰，则平滑版本将在该峰周围下降到原始版本以下。然后，我可以简单地使用此内容和原始内容的最大值，但这会在转换发生时引入非平滑斑点。

此外，我希望能够使用前瞻和后瞻来为此结果曲线参数化算法，因此，如果具有大前瞻和小后瞻，则生成的曲线将更倾向于下降边缘，而具有大后瞻和小前瞻，则将更接近上升边缘。

我尝试使用pandas。 系列（a）.rolling（）设施获得滚动平均值，滚动最大值等，但到目前为止，我找不到一种方法来生成我的输入的平滑版本，该版本在所有情况下都保持在输入之上。

我猜测有一种方法可以将滚动最大值和滚动平均值结合起来以实现我想要的效果，因此这里提供了一些计算这些的代码：

import pandas as pd
import numpy as np

我的输入曲线：

original = np.array([ 5, 5, 5, 8, 8, 8, 2, 2, 2, 2, 2, 3, 3, 7 ])

这可以通过填充左侧（pre）和右侧（post）的边缘值来为任何滚动函数做准备：

pre = 2
post = 3
padded = np.pad(original, (pre, post), 'edge')

现在我们可以应用滚动平均：

smoothed = pd.Series(padded).rolling(
    pre + post + 1).mean().get_values()[pre+post:]

但现在平滑版本在原始版本下方，例如在索引4处：

print(original[4], smoothed[4])  # 8 and 5.5

计算滚动最大值，你可以使用这个：

maximum = pd.Series(padded).rolling(
    pre + post + 1).max().get_values()[pre+post:]

但是仅仅使用滚动最大值在许多情况下并不平滑，并且会在原始峰值周围显示许多平顶。我更喜欢一种平滑处理这些峰值的方法。

如果您还安装了pyqtgraph，您可以轻松地绘制这样的曲线：

import pyqtgraph as pg
p = pg.plot(original)
p.plotItem.plot(smoothed, pen=(255,0,0))

当然，其他绘图库也可以。

我想要的结果是一条曲线，例如由这些值形成的曲线：

goal = np.array([ 5, 7, 7.8, 8, 8, 8, 7, 5, 3.5, 3, 4, 5.5, 6.5, 7 ])

这是曲线的图像。白线是原始数据（输入），红色是滚动平均值，绿色是我想要的大致效果：

编辑：我刚刚发现了一个名为peakutils的模块中的baseline()和envelope()函数。这两个函数可以计算给定度数的多项式，分别适用于输入的下峰值和上峰值。对于小样本来说，这可能是一个不错的解决方案。但如果应用于数百万个值的大样本，则需要非常高的度数，这时计算时间也相当长。将其分段处理（每段一部分）会带来一堆新的问题和难题（如如何正确拼接并保持平滑且保证在输入之上，处理大量碎片时的性能等），所以如果可能的话，我想避免这种情况。

编辑2：我有一个很有前途的方法，即通过反复应用一个过滤器来创建滚动均值，稍微向左和向右移动它，然后取这两个和原始样本的最大值。经过多次应用后，它可以以我想要的方式平滑曲线。然而，在深谷中可能仍然存在一些不平滑的地方。以下是此代码：

pre = 30
post = 30
margin = 10
s = [ np.array(sum([[ x ] * 100 for x in
      [ 5, 5, 5, 8, 8, 8, 2, 2, 2, 2, 2, 3, 3, 7 ]], [])) ]
for _ in range(30):
  s.append(np.max([
    pd.Series(np.pad(s[-1], (margin+pre, post), 'edge')).rolling(
      1 + pre + post).mean().get_values()[pre+post:-margin],
    pd.Series(np.pad(s[-1], (pre, post+margin), 'edge')).rolling(
      1 + pre + post).mean().get_values()[pre+post+margin:],
    s[-1]], 0))

这将创建30次应用过滤器的迭代，使用pyqtplot可以绘制这些迭代的图形：

p = pg.plot(original)
for q in s:
  p.plotItem.plot(q, pen=(255, 100, 100))

生成的图像如下所示：

我不喜欢这种方法的两个方面：①需要重复多次（这会减慢我的速度），②在山谷部分仍然存在不平滑的部分（尽管在我的用例中可能可以接受）。

- Alfe

---

嗯，绿色曲线只是手绘示意图。正如我所写的，它理想上应该可配置为更倾向于左侧或右侧斜坡。我认为手绘的这个似乎在避开上升的斜坡并倾向于下降的斜坡。 - Alfe

以下与上面相同（只是在y轴上反转信号），所以你要找的是同样的东西。我想到的最好的解决方案还没有机会在这里讨论：创建一个滚动最大值，然后使用与窗口宽度相同的高斯窗口来创建一个滚动平均值。结果是一个光滑的曲线，严格位于输入之上，而且完全对称；我希望左右两边有不同的斜率。你考虑的圆可能会陷入狭窄的沟壑并撞上墙壁，导致曲线不平滑。稍后我会详细说明我的方法。 - Alfe

是的，有一点。如果它们是尖锐的（那么它们本来就应该被平滑），这是最明显的，但是已经平滑的峰值也会变得稍微宽一些。我很快会添加代码和图表。 - Alfe

@endolith 请看一下我的答案。它最接近我所寻找的内容。也许那也符合你的需求。 - Alfe

@Alfe 我想我应该单独提一个问题 :) - endolith

Answer 1

我已经玩过了一段时间，我认为我找到了两个解决我直接需求的主要答案。我将在下面给出它们。

import numpy as np
import pandas as pd
from scipy import signal
import pyqtgraph as pg

以下是必要的导入，用于所有代码。当然，pyqtgraph 仅用于显示内容，因此您实际上不需要它。

对称平滑

这可以用于创建总是在信号之上的平滑线，但它不能区分上升和下降的边缘，因此单个峰值周围的曲线将看起来是对称的。在许多情况下，这可能很好，因为它比下面的非对称解决方案简单得多（并且也没有我所知道的任何怪癖）。

s = np.repeat([5, 5, 5, 8, 8, 8, 2, 2, 2, 2, 2, 3, 3, 7], 400) + 0.1
s *= np.random.random(len(s))
pre = post = 400
x = pd.Series(np.pad(s, (pre, post), 'edge')).rolling(
    pre + 1 + post).max().get_values()[pre+post:]
y = pd.Series(np.pad(x, (pre, post), 'edge')).rolling(
    pre + 1 + post, win_type='blackman').mean().get_values()[pre+post:]
p = pg.plot(s, pen=(100,100,100))
for c, pen in ((x, (0, 200, 200)),
               (y, pg.mkPen((255, 255, 255), width=3, style=3))):
    p.plotItem.plot(c, pen=pen)

• 创建一个滚动最大值 (x, 青色), 并且
• 创建一个窗口化的滚动均值 (z, 白色虚线)。

非对称平滑处理

我的使用情况需要一种版本，可以区分上升和下降边缘。输出速度在下降或上升时应该是不同的。

评论: 用作压缩器/扩展器的包络时，快速上升的曲线意味着几乎完全抑制突然的噪声效果，而缓慢上升的曲线意味着在响声出现之前长时间压缩信号，保持动态。另一方面，如果曲线在响声后迅速下降，这将使得响声后的安静部分变得更加清晰可听，而缓慢下降的曲线也会保持动态，并仅缓慢地将信号扩展回正常水平。

s = np.repeat([5, 5, 5, 8, 8, 8, 2, 2, 2, 2, 2, 3, 3, 7], 400) + 0.1
s *= np.random.random(len(s))
pre, post = 100, 1000
t = pd.Series(np.pad(s, (post, pre), 'edge')).rolling(
    pre + 1 + post).max().get_values()[pre+post:]
g = signal.get_window('boxcar', pre*2)[pre:]
g /= g.sum()
u = np.convolve(np.pad(t, (pre, 0), 'edge'), g)[pre:]
g = signal.get_window('boxcar', post*2)[:post]
g /= g.sum()
v = np.convolve(np.pad(t, (0, post), 'edge'), g)[post:]
u, v = u[:len(v)], v[:len(u)]
w = np.min(np.array([ u, v ]),0)
pre = post = max(100, min(pre, post)*3)
x = pd.Series(np.pad(w, (pre, post), 'edge')).rolling(
    pre + 1 + post).max().get_values()[pre+post:]
y = pd.Series(np.pad(x, (pre, post), 'edge')).rolling(
    pre + 1 + post, win_type='blackman').mean().get_values()[pre+post:]
p = pg.plot(s, pen=(100,100,100))
for c, pen in ((t, (200, 0, 0)),
               (u, (200, 200, 0)),
               (v, (0, 200, 0)),
               (w, (200, 0, 200)),
               (x, (0, 200, 200)),
               (y, pg.mkPen((255, 255, 255), width=3))):
    p.plotItem.plot(c, pen=pen)

这个序列无情地结合了几种信号处理方法。

• 输入信号显示为灰色。它是上面提到的输入的嘈杂版本。
• 对此应用滚动最大值 (t，红色)。
• 然后创建一个特殊设计的下降边缘卷积曲线 (g)，并计算卷积 (u，黄色)。
• 使用不同的卷积曲线 (再次使用 g) 重复上述过程以获取上升边缘，并计算卷积 (v，绿色)。
• u 和 v 的最小值是具有所需斜率但仍不太平滑的曲线；特别是当下降和上升斜率相互到达时，它会产生丑陋的尖峰 (w，紫色)。
• 在此基础上应用上述对称方法：
  • 创建此曲线的滚动最大值 (x，青色)。
  • 创建此曲线的窗口化滚动均值 (y，白色虚线)。

- Alfe

Answer 2

正如我在注释中指出的那样，您的绿线在八高台之前和之后的区域中表现不一致。如果您想要左侧区域的行为，可以尝试以下方法：

import numpy as np 
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.spatial import ConvexHull 

# %matplotlib inline # for interactive notebooks

y=np.array([ 5, 5, 5, 8, 8, 8, 2, 2, 2, 2, 2, 3, 3, 7])
x=np.array(range(len(y)))
#######
# This essentially selects the vertices that you'd touch streatching a 
# rubber band over the top of the function
vs = ConvexHull(np.asarray([x,y]).transpose()).vertices 
indices_of_upper_hull_verts = list(reversed(np.concatenate([vs[np.where(vs == len(x)-1)[0][0]: ],vs[0:1]])))

newX = x[indices_of_upper_hull_verts]
newY = y[indices_of_upper_hull_verts]
#########

x_smooth = np.linspace(newX.min(), newX.max(),500)
f = interp1d(newX, newY, kind='quadratic')
y_smooth=f(x_smooth)

plt.plot (x,y)
plt.plot (x_smooth,y_smooth)
plt.scatter (x, y)

这将产生：

更新：

这里有一个替代方案，可能更适合您。如果您使用以1为中心的简单卷积代替滚动平均值，则得到的曲线始终大于输入。卷积核的翼部可以调整为向前/向后查看。

像这样：

import numpy as np 
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.ndimage.filters import convolve

## For interactive notebooks
#%matplotlib inline 

y=np.array([ 5, 5, 5, 8, 8, 8, 2, 2, 2, 2, 2, 3, 3, 7]).astype(np.float)

preLength = 1
postLength = 1
preWeight = 0.2
postWeight = 0.2

kernal = [preWeight/preLength for i in range(preLength)] + [1] + [postWeight/postLength for i in range(postLength)]

output = convolve(y,kernal)

x=np.array(range(len(y)))

plt.plot (x,y)
plt.plot (x,output)

plt.scatter (x, y)

一个缺点是，由于集成的内核通常比1大（这确保了输出曲线平滑且永远不会低于输入），因此输出曲线将始终大于输入曲线，例如在大峰值之上而不是像您所画的那样完全重合。

- Joshua R.

---

我今天稍后会研究一下。乍一看，它似乎过于平滑，根本没有掉入大山谷。但也许当我更仔细地分析答案时，我会发现更多信息。 - Alfe

@Alfe，我添加了一个使用scipy.ndimage.filters.convolve的替代方案-它可能更好地满足您的需求。 - Joshua R.

是的，我忘了提到我想触及原始峰值，而不是达到比必要更高的峰值。将整个输入曲线的最大值作为常量值也会给我一个非常平滑的曲线（常数的导数为0），但当然这不是我想要的...但这是另一种方法，也许我们可以将其与其他东西结合起来 :) - Alfe

Answer 3

作为解决问题的初步尝试，我编写了一个函数，通过最小化积分来将多项式拟合到数据上，并受到多项式严格高于数据点的约束。我猜想如果你在数据上逐段进行这个过程，它可能对你有用。

import scipy.optimize

def upperpoly(xdata, ydata, order):
    def objective(p):
        """Minimize integral"""
        pint = np.polyint(p)
        integral = np.polyval(pint, xdata[-1]) - np.polyval(pint, xdata[0])
        return integral

    def constraints(p):
        """Polynomial values be > data at every point"""
        return np.polyval(p, xdata) - ydata

    p0 = np.polyfit(xdata, ydata, order)
    y0 = np.polyval(p0, xdata)
    shift = (ydata - y0).max()
    p0[-1] += shift

    result = scipy.optimize.minimize(objective, p0, 
                                     constraints={'type':'ineq', 
                                                  'fun': constraints})

    return result.x

- chthonicdaemon

---

我的目标是将其应用于由数百万个值组成的曲线。使用纯Python迭代（为了应用分段算法而进行数千次迭代），而不是一次性在整个样本上使用Numpy、Scipy或Pandas函数，我担心速度会太慢。但是我稍后会仔细研究这个问题。 - Alfe

请查看滚动球基线去除算法。 - chthonicdaemon