作者：CSDN @ _养乐多_

在科学研究中，经常需要对实验数据进行拟合，以找出其中的规律。本文介绍了如何使用Python中的NumPy、Matplotlib和SciPy库进行谐波拟合，并对拟合结果进行残差分析。

从下图可以看出，谐波拟合曲线（红色线）很好地拟合了原始数据（蓝色线），表明所选择的谐波拟合函数能够很好地描述数据的变化趋势。
残差序列（绿色线）在零线附近波动，说明拟合效果较好，残差没有明显的系统性偏差。

---

---

一、谐波拟合

1.1 谐波拟合

谐波拟合是一种拟合方法，用于拟合具有周期性结构的数据。在谐波拟合中，我们假设观测数据可以由一个或多个谐波函数的线性组合来描述。谐波函数是正弦或余弦函数，其周期性特征使得它们适合拟合周期性数据，比如周期性信号或周期性震动。

谐波拟合的公式通常是一个或多个谐波函数的线性组合。对于单个谐波函数，其一般形式可以表示为：

$$f(x)=Asin(2\pi \omega x+\varphi )+C$$

其中：

• A 是振幅（Amplitude）

• ω 是角频率（Angular frequency），与周期 T 之间的关系是 ω=2π​ \ T

• ϕ 是相位（Phase）

• C 是偏移（Offset）

多个谐波函数的线性组合可以通过对单个谐波函数的形式进行叠加得到。

例如，如果我们要拟合两个谐波函数的线性组合，则谐波拟合的公式可以表示为：

$$f(x)=A1sin(2\pi \omega 1x+\varphi 1)+A2sin(2\pi \omega 2x+\varphi 2)+C$$

其中 A1、A2、ω1、ω2、ϕ1、ϕ1、C 分别是两个谐波函数的振幅、角频率、相位和偏移，拟合过程就是通过最小二乘法来确定这些参数的值，使得拟合曲线与实际观测数据的残差的平方和最小化。

1.2 最小二乘法

最小二乘法是一种常用的参数估计方法，用于找到一组参数，使得拟合曲线与观测数据的残差的平方和最小化。在谐波拟合中，最小二乘法用于确定谐波函数的振幅、周期、相位和偏移等参数，以使谐波拟合曲线与观测数据的差异最小化。

1.3 最小二乘法给谐波拟合确定参数的步骤

步骤如下：

1. 定义谐波拟合函数：首先，需要定义一个谐波拟合函数，其形式通常是一个或多个谐波函数的线性组合。这个函数的参数包括谐波函数的振幅、周期、相位和偏移等参数。

2. 构造示例数据：准备需要拟合的实际观测数据，这些数据可能具有周期性结构。

3. 初始化参数估计值：为谐波拟合函数的参数提供一个初始估计值，这些估计值可以根据数据的特点和先验知识来确定。

4. 使用最小二乘法进行拟合：调用最小二乘法的算法，如 scipy.optimize.curve_fit，将谐波拟合函数、实际数据和初始参数估计值作为输入，通过最小化残差的平方和来确定参数的最优值。

5. 获取拟合参数：最小二乘法确定参数后，将得到谐波拟合函数的振幅、周期、相位和偏移等参数。

6. 绘制拟合曲线：使用得到的参数值，将谐波拟合函数应用于原始数据，绘制拟合曲线。

通过这些步骤，最小二乘法可以帮助我们找到最优的参数值，从而得到最佳的谐波拟合结果，使得拟合曲线与实际观测数据的拟合效果最好。

二、完整代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义谐波拟合函数
def harmonic_fit(x, amplitude, period, phase, offset):
    """谐波拟合函数"""
    return amplitude * np.sin(2 * np.pi * x / period + phase) + offset

# 构造示例数据
x_data = np.linspace(0, 10, 100)
y_data = 5 * np.sin(2 * np.pi * x_data / 3 + np.pi / 4) + 2 + np.random.normal(0, 0.5, 100)

# 进行谐波拟合
initial_guess = [1, 3, np.pi / 4, 0]  # 初始参数估计值
params, covariance = curve_fit(harmonic_fit, x_data, y_data, p0=initial_guess)

# 计算谐波拟合曲线和残差序列
fit_curve = harmonic_fit(x_data, *params)
residuals = y_data - fit_curve

# 绘图
plt.figure(figsize=(8, 6))

# 绘制原始数据和谐波拟合曲线
plt.subplot(211)  # 上半部分绘制原始数据和谐波拟合曲线
plt.plot(x_data, y_data, label='原始数据', color='blue')
plt.plot(x_data, fit_curve, label='谐波拟合', color='red')
plt.xlabel('时间')
plt.ylabel('幅值')
plt.title('原始数据与谐波拟合曲线')
plt.legend()

# 绘制残差序列
plt.subplot(212)  # 下半部分绘制残差序列
plt.plot(x_data, residuals, label='残差序列', color='green')
plt.xlabel('时间')
plt.ylabel('残差')
plt.title('残差序列')
plt.legend()

plt.tight_layout()  # 自动调整子图间距
plt.show()