在Python中平滑曲线，同时保留端点处的值和坡度

Question

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

非常感谢你的帮助

附言：所使用的种子是一种具有挑战性的种子

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

Answer 1

是的，最小化是处理此平滑问题的一个好方法

最小二乘问题

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

最小（1/2）和（u[n]-s[n]）²+（λ/2）和（u[n+1]-2u[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]-4u[n+1]+6u[n]-4u[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不是太大，或者如果N很大，可以使用numpy.linalg.solve来执行此操作，或者尝试一种类似共轭梯度的迭代方法

Answer 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

可以为多个平滑度值绘制此曲线，并拾取适合您需要的曲线。通过定义起点和终点，也可以仅对实际曲线中的一系列值应用此方法