在 Python 中平滑曲线，同时保留端点处的值和斜率
Smooth a curve in Python while preserving the value and slope at the end points
其实我有两个解决这个问题的方法，它们都在下面应用到一个测试用例中。问题是 none 是完美的：第一个只考虑了两个端点，另一个不能“任意平滑”：一个人可以达到的平滑度是有限的（我正在展示的那个）。
我确信有更好的解决方案，从第一个解决方案到另一个解决方案，一直到完全不平滑。它可能已经在某处实施。也许用任意数量的样条等分布解决最小化问题？
非常感谢您的帮助
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()
是的，最小化是解决这个平滑问题的好方法。
最小二乘问题
这里有一个最小二乘公式的建议：令 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] (斜率约束),
其中在最小化 objective 中，总和超过 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，将已知量移动到 right-hand 侧：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 很大，请尝试使用共轭梯度之类的迭代方法。
您可以应用一种简单的平滑方法并绘制具有不同平滑度值的平滑曲线，以查看哪种效果最好。
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
您可以为多个平滑度值绘制此曲线，并选择适合您需要的曲线。您还可以通过定义 start 和 end
将此方法仅应用于实际曲线中的一系列值
其实我有两个解决这个问题的方法，它们都在下面应用到一个测试用例中。问题是 none 是完美的：第一个只考虑了两个端点，另一个不能“任意平滑”：一个人可以达到的平滑度是有限的（我正在展示的那个）。 我确信有更好的解决方案，从第一个解决方案到另一个解决方案，一直到完全不平滑。它可能已经在某处实施。也许用任意数量的样条等分布解决最小化问题？
非常感谢您的帮助
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()
是的，最小化是解决这个平滑问题的好方法。
最小二乘问题
这里有一个最小二乘公式的建议：令 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] (斜率约束),
其中在最小化 objective 中，总和超过 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，并对
u[2:-2] -= (0.05 / λ) * (u - s)[2:-2] + np.convolve(u, [1, -4, 6, -4, 1])[4:-4]
但是通过更多的工作，可以通过直接求解 E-L 方程来做得更好。对于每个 n，将已知量移动到 right-hand 侧：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 很大，请尝试使用共轭梯度之类的迭代方法。
您可以应用一种简单的平滑方法并绘制具有不同平滑度值的平滑曲线，以查看哪种效果最好。
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
您可以为多个平滑度值绘制此曲线，并选择适合您需要的曲线。您还可以通过定义 start 和 end
将此方法仅应用于实际曲线中的一系列值