用Python和NumPy手把手教你：从方波合成动画看懂傅里叶级数（附完整代码）

用Python和NumPy手把手教你：从方波合成动画看懂傅里叶级数（附完整代码）

傅里叶级数这个数学概念，听起来总是带着几分神秘色彩。但当我第一次看到用动画展示方波如何由不同频率的正弦波叠加而成时，那种直观的感受彻底改变了我对这个抽象概念的理解。本文将带你用Python亲手实现这个神奇的合成过程，通过代码和可视化，让傅里叶级数从数学公式变成生动的视觉体验。

1. 环境准备与基础概念

在开始编码之前，我们需要确保环境配置正确。这个项目只需要三个Python库：

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation

傅里叶级数的核心思想是：任何周期函数都可以表示为不同频率正弦和余弦函数的无限和。对于周期为2π的方波函数，其数学表达式为：

def square_wave(t): return np.where((t % (2*np.pi)) < np.pi, 1, -1)

这个函数会在0到π区间输出1，在π到2π区间输出-1，如此循环往复，形成一个标准的方波。

傅里叶级数告诉我们，这个方波可以分解为：

4/π * [sin(t) + (1/3)sin(3t) + (1/5)sin(5t) + ...]

这个级数只包含奇数次的正弦波，每个分量的振幅随着频率增加而减小。理解这一点对后续的代码实现至关重要。

2. 构建傅里叶级数近似函数

让我们从实现傅里叶级数的部分和开始。定义一个函数来计算包含n项谐波的方波近似：

def fourier_series(t, n_terms): result = np.zeros_like(t) for n in range(1, n_terms*2, 2): # 只考虑奇数项 result += (4/(np.pi*n)) * np.sin(n*t) return result

为了理解这个近似过程的效果，我们可以绘制不同谐波数量下的近似曲线：

t = np.linspace(0, 4*np.pi, 1000) plt.figure(figsize=(10,6)) for n in [1, 3, 5, 10, 20, 50]: plt.plot(t, fourier_series(t, n), label=f'{n} terms') plt.plot(t, square_wave(t), 'k--', label='True square wave') plt.legend() plt.show()

随着谐波数量的增加，近似曲线会越来越接近理想的方波，但在跳变点附近会出现所谓的"吉布斯现象"——过冲不会消失，只是变窄。

3. 创建动态合成动画

静态图像难以展示傅里叶级数的动态合成过程，因此我们创建一个动画来直观展示每个正弦波如何贡献到最终的方波形状：

def create_animation(): fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,8)) t = np.linspace(0, 4*np.pi, 1000) line, = ax1.plot(t, np.zeros_like(t), 'r') ax1.plot(t, square_wave(t), 'k--') ax1.set_ylim(-1.5, 1.5) # 初始化所有谐波 harmonics = [] for n in range(1, 20, 2): harmonic, = ax2.plot(t, (4/(np.pi*n))*np.sin(n*t), alpha=0.5) harmonics.append(harmonic) def update(frame): n_terms = frame + 1 # 更新合成波 line.set_ydata(fourier_series(t, n_terms)) # 高亮当前添加的谐波 for i, h in enumerate(harmonics): h.set_alpha(1.0 if i == frame else 0.3) ax1.set_title(f'Square wave approximation with {n_terms} terms') return [line] + harmonics ani = FuncAnimation(fig, update, frames=10, interval=1000, blit=True) plt.tight_layout() return ani

这个动画会逐步添加每个谐波分量，同时在上方面板显示合成结果，在下方面板高亮当前添加的谐波。你可以保存为GIF或直接在Jupyter中显示：

ani = create_animation() ani.save('fourier_animation.gif', writer='pillow', fps=1)

4. 深入理解吉布斯现象

当我们在代码中增加谐波数量时，会发现一个有趣的现象：在方波的跳变点附近，合成波形会出现过冲和振荡，即使增加更多谐波，这个过冲也不会完全消失，只是变得更加集中在跳变点附近。这就是著名的吉布斯现象。

我们可以量化分析这个现象：

def analyze_gibbs(n_terms_list): t = np.linspace(0, 2*np.pi, 10000) jump_point = np.pi idx = np.argmin(np.abs(t - jump_point)) overshoots = [] for n in n_terms_list: approx = fourier_series(t, n) max_val = np.max(approx[idx-50:idx+50]) overshoots.append((max_val - 1) * 100) # 百分比过冲 plt.plot(n_terms_list, overshoots, 'o-') plt.xlabel('Number of terms') plt.ylabel('Overshoot percentage') plt.title('Gibbs phenomenon analysis') plt.show() analyze_gibbs(range(1, 101, 5))

数学上可以证明，这个过冲大约为跳变高度的9%，无论添加多少谐波都不会消除。这是傅里叶级数在间断点收敛的一个特性。

5. 性能优化与交互式探索

当我们需要计算大量谐波时，原始的实现可能会变得缓慢。我们可以利用NumPy的向量化运算来优化：

def optimized_fourier(t, n_terms): n = np.arange(1, 2*n_terms, 2) # 所有奇数 coeffs = 4 / (np.pi * n) # 对应系数 sins = np.sin(np.outer(t, n)) # 所有正弦项 return np.dot(sins, coeffs) # 向量化计算

为了更灵活地探索不同参数的影响，我们可以创建一个交互式工具：

from ipywidgets import interact def interactive_exploration(n_terms=(1, 50, 2)): t = np.linspace(0, 4*np.pi, 1000) plt.figure(figsize=(10,5)) plt.plot(t, square_wave(t), 'k--', label='True square wave') plt.plot(t, optimized_fourier(t, n_terms), 'r', label=f'{n_terms} terms') plt.ylim(-1.5, 1.5) plt.legend() plt.show() interact(interactive_exploration)

这个交互式工具让你可以滑动调整谐波数量，实时观察合成效果的变化，非常适合教学和自学使用。

6. 扩展到其他波形

理解了方波的合成后，我们可以将同样的方法应用到其他周期波形上。例如，三角波的傅里叶级数为：

def triangle_wave(t): return (2/np.pi) * np.arcsin(np.sin(t)) def fourier_triangle(t, n_terms): result = np.zeros_like(t) for n in range(1, n_terms*2, 2): result += ((-1)**((n-1)/2)) * (8/(np.pi**2 * n**2)) * np.sin(n*t) return result

锯齿波的傅里叶级数则是：

def fourier_sawtooth(t, n_terms): result = np.zeros_like(t) for n in range(1, n_terms+1): result += ((-1)**(n+1)) * (2/(np.pi*n)) * np.sin(n*t) return result

通过修改我们的动画代码，可以轻松创建这些波形的合成动画，比较它们收敛特性的差异。

7. 实际应用与高级话题

傅里叶级数不仅仅是数学上的奇观，它在信号处理、音频合成、图像压缩等领域有广泛应用。例如，在音频合成中：

def play_wave(freq, duration=2, sample_rate=44100, wave_type='sine'): t = np.linspace(0, duration, int(sample_rate*duration), False) if wave_type == 'sine': signal = np.sin(2*np.pi*freq*t) elif wave_type == 'square': signal = np.sign(np.sin(2*np.pi*freq*t)) elif wave_type == 'triangle': signal = (2/np.pi) * np.arcsin(np.sin(2*np.pi*freq*t)) # 播放或保存音频...

理解傅里叶级数还能帮助我们更好地掌握傅里叶变换，这是信号频域分析的基础。在NumPy中，我们可以用np.fft模块进行快速傅里叶变换分析。

通过这个项目，我们不仅实现了方波的傅里叶合成，还建立了一套可以扩展到其他波形分析的工具。这种从数学到代码再到可视化的完整过程，是理解复杂概念的强大方法。