移动最小二乘法拟合趋势面滤波代码

以下是一个基于移动最小二乘法的趋势面滤波的代码示例：

import numpy as np
from numpy.linalg import inv

def moving_least_squares(data, window_size):
    """
    移动最小二乘法拟合趋势面滤波
    :param data: 输入数据，形状为 (n, 2)，其中第一列是横坐标，第二列是纵坐标
    :param window_size: 窗口大小，必须为奇数
    :return: 滤波后的数据，形状与输入数据相同
    """
    assert window_size % 2 == 1, "窗口大小必须为奇数"

    n = data.shape[0]
    result = np.zeros_like(data)

    for i in range(n):
        # 计算当前窗口的索引
        indices = np.arange(max(0, i - window_size // 2), min(n, i + window_size // 2 + 1))
        m = len(indices)

        # 构造矩阵 A 和向量 b
        A = np.zeros((m, 3))
        b = np.zeros((m, 1))

        for j in range(m):
            x, y = data[indices[j], :]
            A[j, 0] = x ** 2
            A[j, 1] = x
            A[j, 2] = 1
            b[j, 0] = y

        # 计算最小二乘解
        x = inv(A.T @ A) @ A.T @ b

        # 将拟合结果保存到结果数组中
        result[i, 0] = data[i, 0]
        result[i, 1] = A[m // 2, 0] * x[0, 0] + A[m // 2, 1] * x[1, 0] + A[m // 2, 2] * x[2, 0]

    return result

使用方法如下：

import matplotlib.pyplot as plt

# 生成测试数据
x = np.linspace(0, 2 * np.pi, 100)
y = np.sin(x) + 0.1 * np.random.randn(100)

# 拼接成二维数组
data = np.vstack([x, y]).T

# 执行移动最小二乘法滤波
filtered_data = moving_least_squares(data, 11)

# 绘制原始数据和滤波结果
plt.plot(x, y, label='Original')
plt.plot(filtered_data[:, 0], filtered_data[:, 1], label='Filtered')
plt.legend()
plt.show()

可以看到，使用移动最小二乘法可以有效地去除原始数据中的噪声，还原出数据的趋势面。