Python中三角脉冲的傅里叶级数拟合问题

问题根源分析 #

原代码的核心错误在于傅里叶级数重构时的频率项计算错误，以及FFT系数的归一化缺失：

• 未考虑信号的周期T=20，基频应为$\omega_0 = 2\pi/T = \pi/10$，重构时的指数项应该是$\exp(1j \cdot n \cdot \omega_0 \cdot t)$，而非直接用$n \cdot t$。
• np.fft.fft输出的系数未做归一化，需要除以FFT点数N才能得到正确的傅里叶级数系数。
• 手动遍历所有FFT系数时，包含了正、负频率分量（FFT的对称性），直接累加会导致重复计算，需合理处理或使用逆FFT简化流程。

修正方案 #

方案1：使用np.fft.ifft直接重构（最简单可靠） #

利用numpy内置的逆FFT函数，自动处理频率和归一化问题：

import numpy as np
import matplotlib.pyplot as plt

def Real(t):
    if t < 0:
        return 298
    elif t < 10:
        return 10.2*t + 298
    elif t < 20:
        return 502 - 10.2*t
    else:
        return 298

# 生成原信号数据
T = 20  # 信号周期
ts = np.linspace(0, T, 100, endpoint=False)  # 避免重复周期端点
Rea = np.array([Real(t) for t in ts])

N = 200  # FFT点数
Coeff = np.fft.fft(Rea, n=N)
# 逆FFT重构，自动处理归一化和频率
Fourier = np.fft.ifft(Coeff)
# 取实部并对齐时间轴（补零后的FFT对应更长的时间序列，取前100点匹配原信号）
Fourier = np.real(Fourier[:len(ts)])

# 绘图对比
plt.plot(ts, Rea, label='原信号')
plt.plot(ts, Fourier, label='FFT重构信号', linestyle='--')
plt.legend()
plt.show()

方案2：手动实现傅里叶级数重构（理解原理） #

手动处理频率项和归一化，完整遍历所有FFT分量：

import numpy as np
import matplotlib.pyplot as plt

def Real(t):
    if t < 0:
        return 298
    elif t < 10:
        return 10.2*t + 298
    elif t < 20:
        return 502 - 10.2*t
    else:
        return 298

T = 20
omega0 = 2 * np.pi / T  # 基频
ts = np.linspace(0, T, 100, endpoint=False)
Rea = np.array([Real(t) for t in ts])

N = 200
Coeff = np.fft.fft(Rea, n=N) / N  # 归一化系数

def Four_Func(t):
    ans = 0
    for n in range(N):
        # 正确的频率项：n*omega0*t
        ans += Coeff[n] * np.exp(1j * n * omega0 * t)
    return np.real(ans)

Fourier = np.array([Four_Func(t) for t in ts])

# 绘图对比
plt.plot(ts, Rea, label='原信号')
plt.plot(ts, Fourier, label='手动重构信号', linestyle='--')
plt.legend()
plt.show()

额外优化点 #

• 生成时间序列时使用endpoint=False，避免周期端点重复采样，减少FFT的边界误差。
• 如果需要仅保留正弦和余弦项（而非复指数），可以将复系数转换为实部（余弦项系数）和虚部（正弦项系数）：

  # 直流分量
  a0 = 2 * np.real(Coeff[0])
  # 余弦项系数a_n，正弦项系数b_n
  a_n = 2 * np.real(Coeff[1:N//2])
  b_n = -2 * np.imag(Coeff[1:N//2])
  
  def Four_Real_Func(t):
      ans = a0 / 2
      for n in range(1, N//2):
          ans += a_n[n-1] * np.cos(n*omega0*t) + b_n[n-1] * np.sin(n*omega0*t)
      return ans
  

内容的提问来源于stack exchange，提问作者KNVCSG