在 Python 中平滑曲线，同时保留端点的值和斜率

Posted 2023-03-10

【中文标题】在 Python 中平滑曲线，同时保留端点的值和斜率

【英文标题】：Smooth a curve in Python while preserving the value and slope at the end points

【发布时间】：2021-02-04 20:09:29

【问题描述】：

我实际上有两个解决这个问题的方法，它们都适用于下面的测试用例。问题是它们都不是完美的：第一个只考虑两个端点，另一个不能“任意平滑”：一个可以达到的平滑量是有限的（我正在显示）。 我确信有一个更好的解决方案，从第一个解决方案到另一个解决方案，一直到根本没有平滑。它可能已经在某处实施。也许用任意数量的样条线等分布来解决最小化问题？

非常感谢您的帮助

Ps：使用的种子很有挑战性

    import matplotlib.pyplot as plt
    from scipy import interpolate
    from scipy.signal import savgol_filter
    import numpy as np 
    import random

    def scipy_bspline(cv, n=100, degree=3):
        """ Calculate n samples on a bspline
            cv :      Array ov control vertices
            n  :      Number of samples to return
            degree:   Curve degree
        """
        cv = np.asarray(cv)
        count = cv.shape[0]
        degree = np.clip(degree,1,count-1)
        kv = np.clip(np.arange(count+degree+1)-degree,0,count-degree)
    
        # Return samples
        max_param = count - (degree * (1-periodic))
        spl = interpolate.BSpline(kv, cv, degree)
        return spl(np.linspace(0,max_param,n))

    def round_up_to_odd(f):
        return np.int(np.ceil(f / 2.) * 2 + 1)
    
    def generateRandomSignal(n=1000, seed=None):
        """
        Parameters
        ----------
        n : integer, optional
            Number of points in the signal. The default is 1000.
    
        Returns
        -------
        sig : numpy array
    
        """
        np.random.seed(seed)
        print("Seed was:", seed)
        steps = np.random.choice(a=[-1, 0, 1], size=(n-1))
        roughSig = np.concatenate([np.array([0]), steps]).cumsum(0)
        sig = savgol_filter(roughSig, round_up_to_odd(n/10), 6)
        return sig
    
    # Generate a random signal to illustrate my point
    n = 1000
    t = np.linspace(0, 10, n)
    seed = 45136. # Challenging seed
    sig = generateRandomSignal(n=1000, seed=seed)
    sigInit = np.copy(sig)
    
    # Add noise to the signal
    mean = 0
    std = sig.max()/3.0
    num_samples = n/5
    idxMin = n/2-100
    idxMax = idxMin + num_samples
    tCut = t[idxMin+1:idxMax]
    noise = np.random.normal(mean, std, size=num_samples-1) + 2*std*np.sin(2.0*np.pi*tCut/0.4)
    sig[idxMin+1:idxMax] += noise
    
    # Define filtering range enclosing the noisy area of the signal
    idxMin -= 20
    idxMax += 20
    
    # Extreme filtering solution
    # Spline between first and last points, the points in between have no influence
    sigTrim = np.delete(sig, np.arange(idxMin,idxMax))
    tTrim = np.delete(t, np.arange(idxMin,idxMax))
    f = interpolate.interp1d(tTrim, sigTrim, kind='quadratic')
    sigSmooth1 = f(t)
    
    # My attempt. Not bad but not perfect because there is a limit in the maximum
    # amount of smoothing we can add (degree=len(tSlice) is the maximum)
    # If I could do degree=10*len(tSlice) and converging to the first solution
    # I would be done!
    sigSlice = sig[idxMin:idxMax]
    tSlice = t[idxMin:idxMax]
    cv = np.stack((tSlice, sigSlice)).T
    p = scipy_bspline(cv, n=len(tSlice), degree=len(tSlice))
    tSlice = p.T[0]
    sigSliceSmooth = p.T[1]
    sigSmooth2 = np.copy(sig)
    sigSmooth2[idxMin:idxMax] = sigSliceSmooth
    
    # Plot
    plt.figure()
    plt.plot(t, sig, label="Signal")
    plt.plot(t, sigSmooth1, label="Solution 1")
    plt.plot(t, sigSmooth2, label="Solution 2")
    plt.plot(t[idxMin:idxMax], sigInit[idxMin:idxMax], label="What I'd want (kind of, smoother will be even better actually)")
    plt.plot([t[idxMin],t[idxMax]], [sig[idxMin],sig[idxMax]],"o")
    plt.legend()
    plt.show()
    sys.exit()

【参考方案1】：

是的，最小化是解决此平滑问题的好方法。

最小二乘问题

这里是一个最小二乘公式的建议：让 s[0], ..., s[N] 表示给定信号的 N+1 个样本要平滑，并让 L 和 R 是所需的斜率保留在左右端点。找到平滑后的信号 u[0], ..., u[N] 作为

的最小值

min_u (1/2) sum_n (u[n] - s[n])² + (λ/2) sum_n (u[n+1] - 2 u[n] + u[n-1]) ²

受制于 s[0] = u[0], s[N] = u[N] (值约束), L = u[1] - u[0]，R = u[N] - u[N-1]（斜率约束），

在最小化目标中，总和超过 n = 1，...，N-1，λ 是控制平滑强度的正参数。第一项试图使解接近原始信号，第二项惩罚 u 弯曲以鼓励平滑解。

坡度约束要求 u[1] = L + u[0] = L + s[0] 和 u[N-1] = u[N] - R = s[N] - R. 所以我们可以认为最小化只超过内部样本 u[2], ..., u[N-2].

找到最小化器

最小化器满足欧拉-拉格朗日方程

(u[n] - s[n]) / λ + (u[n+2] - 4 u[n+1] + 6 u[n] - 4 u[n-1] + u[n -2]) = 0 对于 n = 2, ..., N-2。

找到近似解的一种简单方法是通过梯度下降：初始化u = np.copy(s)，设置 u[1] = L + s[0] 和 u[N-1] = s[N] - R，然后执行100 次左右的迭代

u[2:-2] -= (0.05 / λ) * (u - s)[2:-2] + np.convolve(u, [1, -4, 6, -4, 1])[4:-4]

但如果做更多的工作，直接求解 E-L 方程可能会比这做得更好。对于每个 n，将已知量移到右侧：s[n] 以及端点 u[0] = s[0], u[1] = L + s[0], u[N-1 ] = s[N] - R，u[N] = s[N]。你将有一个线性系统“A u = b”，矩阵 A 有类似的行

0, ..., 0, 1, -4, (6 + 1/λ), -4, 1, 0, ..., 0.

最后，求解线性系统以找到平滑信号 u。如果 N 不太大，您可以使用numpy.linalg.solve 来执行此操作，或者如果 N 很大，请尝试使用共轭梯度等迭代方法。

【讨论】：

非常感谢您提供这个非常好的答案！我正在努力。

【参考方案2】：

您可以应用一种简单的平滑方法，并绘制具有不同平滑度值的平滑曲线，看看哪个效果最好。

def smoothing(data, smoothness=0.5):
    last = data[0]
    new_data = [data[0]]
    for datum in data[1:]:
        new_value = smoothness * last + (1 - smoothness) * datum
        new_data.append(new_value)
        last = datum
    return new_data

您可以绘制这条曲线以获得多个平滑度值，然后选择适合您需要的曲线。您也可以通过定义起点和终点将此方法仅应用于实际曲线中的一系列值

【讨论】：

非常感谢您的回答！我只是测试了它，它似乎并没有解决问题。平滑度 = 0 -> 两条曲线相等，平滑度 = 1 -> 两条曲线几乎相等（实际上有点偏移）。也许我错过了什么？

en.wikipedia.org/wiki/Exponential_smoothing 我想这可能会对你有所帮助...