曲线平滑 [ Soymilk ]

曲线平滑

categories
>
others

基础平滑算法
拟合类平滑算法
滤波器类
几何/图形平滑

生成噪声曲线并显示

python

import numpy as np
import matplotlib.pyplot as plt

# ========== 生成并写入带噪声曲线数据 ==========
np.random.seed(0)
x = np.linspace(0, 4 * np.pi, 500)
y_true = np.sin(x)
noise = np.random.normal(0, 0.3, x.shape)
y_noisy = y_true + noise

# 保存数据到文件
np.savetxt("data/noisy_data.csv", np.column_stack((x, y_noisy)), delimiter=",", header="x,y_noisy", comments='')

# ========== 从文件读取数据 ==========
data = np.loadtxt("data/noisy_data.csv", delimiter=",", skiprows=1)
x = data[:, 0]
y_noisy = data[:, 1]

plt.plot(x, y_noisy, label='Noisy', color='gray')
plt.title("Original Noisy Curve")
plt.show()

original noisy

基础平滑算法

算法名	简介	优点	缺点
移动平均（Moving Average）	用相邻点平均代替当前点	简单易实现	会模糊边缘，无法处理非线性数据
加权移动平均（Weighted Moving Average）	给中间点更大权重	平滑更自然	参数依赖经验
中值滤波（Median Filter）	用邻域的中位数代替当前点	抗脉冲噪声强	对高频信号损失较大
高斯滤波（Gaussian Filter）	使用高斯核进行平滑	控制光滑程度好	需要选择合适的 σ

python

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
from scipy.signal import medfilt


data = np.loadtxt("data/noisy_data.csv", delimiter=",", skiprows=1)
x = data[:, 0]
y = data[:, 1]

# 1. 简单移动平均
def moving_average(y, window=5):
    return np.convolve(y, np.ones(window) / window, mode='same')

# 2. 加权移动平均
def weighted_moving_average(y, weights=None):
    if weights is None:
        weights = np.array([1, 2, 3, 2, 1])
    weights = weights / weights.sum()
    return np.convolve(y, weights, mode='same')

# 3. 中值滤波
def median_filter(y, kernel_size=5):
    return medfilt(y, kernel_size=kernel_size)

# 4. 高斯滤波
def gaussian_filter(y, sigma=2):
    return gaussian_filter1d(y, sigma=sigma)

# 应用平滑方法
y_ma = moving_average(y)
y_wma = weighted_moving_average(y)
y_med = median_filter(y)
y_gauss = gaussian_filter(y)

# 画图
plt.figure(figsize=(12, 8))
plt.plot(x, y, label='original', alpha=0.4, linestyle='--')
plt.plot(x, y_ma, label='moving average')
plt.plot(x, y_wma, label='weighted moving average')
plt.plot(x, y_med, label='median filter')
plt.plot(x, y_gauss, label='gaussian filter')

plt.legend()
plt.title('basic smooth')
plt.grid(True)
plt.show()

basic smooth

拟合类平滑算法

算法名	简介	优点	缺点
多项式拟合（Polynomial Fitting）	用曲线整体拟合数据	保留趋势	高阶拟合容易震荡
样条插值（Spline Interpolation）	拟合连续光滑曲线（如 B 样条、三次样条）	平滑度高，常用于 CAD、路径规划	计算复杂度略高
Loess / LOWESS	局部加权回归	抗噪声能力强，适合非线性	参数选择敏感，计算较慢

python

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
import statsmodels.api as sm

# 读取数据
data = np.loadtxt("data/noisy_data.csv", delimiter=",", skiprows=1)
x = data[:, 0]
y = data[:, 1]

# 1. 多项式拟合（Polynomial Fit）
def poly_fit(x, y, degree=4):
coeffs = np.polyfit(x, y, deg=degree)
poly_func = np.poly1d(coeffs)
return poly_func(x)

# 2. 样条插值（Cubic Spline）
def spline_fit(x, y):
cs = CubicSpline(x, y)
return cs(x)

# 3. Loess（局部加权回归）
def loess_fit(x, y, frac=0.2):
lowess = sm.nonparametric.lowess
y_loess = lowess(y, x, frac=frac, return_sorted=False)
return y_loess

# 平滑处理
y_poly = poly_fit(x, y, degree=4)
y_spline = spline_fit(x, y)
y_loess = loess_fit(x, y, frac=0.15)

# 绘图
plt.figure(figsize=(12, 8))
plt.plot(x, y, label='original', linestyle='--', alpha=0.5)
plt.plot(x, y_poly, label='poly')
plt.plot(x, y_spline, label='spline')
plt.plot(x, y_loess, label='loess')

plt.title('fitting smooth')
plt.legend()
plt.grid(True)
plt.savefig('fitted_curve.png') # 保存图像
plt.show()

fitting smooth

滤波器类

算法名	简介	优点	缺点
Savitzky–Golay 滤波器	在滑动窗口内做局部多项式拟合	平滑但保留峰值、斜率	不适合强噪声数据
卡尔曼滤波（Kalman Filter）	适合动态系统中的数据平滑	实时性能好	参数调优复杂
Butterworth / 低通滤波器	信号频域滤波法	控制频率响应	需频域变换（如 FFT）

python

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter, butter, filtfilt
from filterpy.kalman import KalmanFilter

# 读取数据
data = np.loadtxt("data/noisy_data.csv", delimiter=",", skiprows=1)
x = data[:, 0]
y = data[:, 1]

def apply_savgol(y, window_length=11, polyorder=3):
    return savgol_filter(y, window_length, polyorder)

def apply_kalman(y):
    kf = KalmanFilter(dim_x=2, dim_z=1)
    kf.x = np.array([[0.], [0.]])
    kf.F = np.array([[1., 1.], [0., 1.]])
    kf.H = np.array([[1., 0.]])
    kf.P *= 1000.
    kf.R = 5
    kf.Q = np.array([[1., 0.], [0., 1.]]) * 0.01

    estimates = []
    for z in y:
        kf.predict()
        kf.update(z)
        estimates.append(kf.x[0, 0])
    return np.array(estimates)

def apply_butterworth(y, cutoff=0.05, fs=1.0, order=5):
    b, a = butter(order, cutoff, fs=fs, btype='low')
    return filtfilt(b, a, y)

y_sg = apply_savgol(y)
y_kf = apply_kalman(y)
y_bf = apply_butterworth(y)

plt.figure(figsize=(12, 8))
plt.plot(x, y, label="original", linestyle='--', alpha=0.6)
plt.plot(x, y_sg, label="Savitzky–Golay")
plt.plot(x, y_kf, label="Kalman")
plt.plot(x, y_bf, label="Butterworth")

plt.title("filter smooth")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True)
plt.savefig("filter_compare.png")
plt.show()

filter smooth

几何/图形平滑

算法名	应用领域	简介
Chaikin 算法	2D 曲线	迭代方式缩短控制点构造圆滑曲线
Bézier 曲线	图形/动画	精确建模平滑曲线
Catmull-Rom 样条	动画插值	过所有控制点且自然平滑

python

import numpy as np
import matplotlib.pyplot as plt


# Chaikin 曲线细分算法
def chaikin_curve(points, iterations=3):
    for _ in range(iterations):
        new_points = []
        for i in range(len(points) - 1):
            p0 = points[i]
            p1 = points[i + 1]
            q = 0.75 * p0 + 0.25 * p1
            r = 0.25 * p0 + 0.75 * p1
            new_points.extend([q, r])
        points = np.array(new_points)
    return points


# Bézier 曲线（de Casteljau 算法）
def bezier_curve(points, num=100):
    def de_casteljau(pts, t):
        while len(pts) > 1:
            pts = [(1 - t) * pts[i] + t * pts[i + 1] for i in range(len(pts) - 1)]
        return pts[0]

    t_values = np.linspace(0, 1, num)
    curve = np.array([de_casteljau(points.copy(), t) for t in t_values])
    return curve


# Catmull-Rom 样条插值
def catmull_rom_spline(points, num=100):
    def catmull_rom(p0, p1, p2, p3, t):
        t2 = t * t
        t3 = t2 * t
        return 0.5 * ((2 * p1) +
                      (-p0 + p2) * t +
                      (2*p0 - 5*p1 + 4*p2 - p3) * t2 +
                      (-p0 + 3*p1 - 3*p2 + p3) * t3)

    curve = []
    for i in range(1, len(points) - 2):
        for t in np.linspace(0, 1, num):
            pt = catmull_rom(points[i - 1], points[i], points[i + 1], points[i + 2], t)
            curve.append(pt)
    return np.array(curve)


# 读取 data.txt 中的二维点
data = np.loadtxt('data/noisy_data.csv', delimiter=',', skiprows=1)
x, y = data[:, 0], data[:, 1]
points = np.column_stack((x, y))

# 应用三种算法
chaikin_pts = chaikin_curve(points, iterations=3)
bezier_pts = bezier_curve(points)
catmull_pts = catmull_rom_spline(points) if len(points) >= 4 else None

# 绘图
plt.figure(figsize=(12, 6))
plt.plot(points[:, 0], points[:, 1], 'o--', label='Original', linewidth=1)

plt.plot(chaikin_pts[:, 0], chaikin_pts[:, 1], label='Chaikin', linewidth=2)
plt.plot(bezier_pts[:, 0], bezier_pts[:, 1], label='Bezier', linewidth=2)

if catmull_pts is not None:
    plt.plot(catmull_pts[:, 0], catmull_pts[:, 1], label='Catmull-Rom', linewidth=2)

plt.legend()
plt.title('Geometry/Animation Smoothing')
plt.grid(True)
plt.tight_layout()
plt.show()

geometry smooth

others

comment:

Valine
LiveRe
ChangYan