别再死记硬背公式了！用Python+SymPy手把手推导方波傅里叶级数（附完整代码）

用Python+SymPy实战推导方波傅里叶级数：从数学公式到可视化分析

在电力电子和信号处理领域，方波是最基础的波形之一。传统教学中，我们往往需要死记硬背各种傅里叶级数公式，却很少有机会真正理解这些系数是如何计算出来的。本文将带你用Python的SymPy库，从零开始推导方波的傅里叶级数，并通过代码实现可视化分析，让你不仅知其然，更知其所以然。

1. 环境准备与工具介绍

1.1 为什么选择SymPy

SymPy是一个纯Python编写的符号计算库，它能够像Mathematica或Maple一样进行符号运算，但又保持了Python的简洁性。相比数值计算库如NumPy，SymPy可以直接处理数学表达式，保留π、√2等符号形式，这正是推导傅里叶级数所需的特性。

安装SymPy非常简单：

pip install sympy matplotlib numpy

1.2 基础数学概念回顾

傅里叶级数将周期函数表示为正弦和余弦函数的无限和：

f(x) = a₀/2 + Σ(aₙcos(nx) + bₙsin(nx))

其中系数计算公式为：

aₙ = (1/π)∫f(x)cos(nx)dx bₙ = (1/π)∫f(x)sin(nx)dx

提示：对于奇函数，所有aₙ系数为零；对于偶函数，所有bₙ系数为零。选择合适的积分区间可以简化计算。

2. 定义方波函数与对称性分析

2.1 构建符号化方波

我们先定义一个周期为2π的方波函数。在SymPy中，可以使用Piecewise函数表示分段函数：

from sympy import * x, n = symbols('x n', real=True) V = symbols('V', positive=True) # 方波幅值 # 定义周期为2π的方波 square_wave = Piecewise( (V, (x >= 0) & (x < pi)), (-V, (x >= pi) & (x < 2*pi)), (square_wave.subs(x, x - 2*pi), True) # 周期性延拓 )

2.2 奇函数特性验证

我们特意将方波的跳变点设在x=0处，使其成为奇函数。验证奇函数性质：

# 验证f(-x) = -f(x) simplify(square_wave.subs(x, -x) + square_wave) == 0 # 应返回True

由于是奇函数，我们立即可以得出两个结论：

1. 直流分量a₀ = 0
2. 所有余弦系数aₙ = 0

这已经将问题简化了一半！

3. 计算傅里叶系数

3.1 正弦系数bₙ的符号计算

现在只需要计算bₙ系数。根据公式：

b_n = (1/pi) * integrate(square_wave * sin(n*x), (x, 0, 2*pi))

SymPy会自动计算这个积分，结果可能看起来比较复杂。我们可以分步计算：

# 分步计算积分 integral_part1 = integrate(V * sin(n*x), (x, 0, pi)) integral_part2 = integrate(-V * sin(n*x), (x, pi, 2*pi)) b_n = simplify((1/pi) * (integral_part1 + integral_part2))

得到的bₙ表达式为：

bₙ = (V*(2 - 2*cos(πn)))/(πn)

3.2 系数简化与奇偶分析

观察这个结果，我们可以进一步简化：

# 利用cos(πn)的特性简化 b_n_simplified = b_n.subs(cos(pi*n), (-1)**n)

这利用了cos(πn) = (-1)ⁿ的数学性质。现在表达式变为：

bₙ = (2V(1 - (-1)ⁿ))/(πn)

当n为偶数时，1 - (-1)ⁿ = 0；当n为奇数时，1 - (-1)ⁿ = 2。因此：

# 最终简化结果 b_n_final = Piecewise( (0, Eq(n%2, 0)), (4*V/(pi*n), True) )

4. 谐波合成与可视化

4.1 前N项和的计算

现在我们可以构建傅里叶级数的前N项和：

def fourier_series(N): series = 0 for k in range(1, N+1): n = 2*k - 1 # 只取奇数 series += b_n_final.subs(n, 2*k-1) * sin((2*k-1)*x) return series

4.2 使用Matplotlib可视化

让我们绘制前1、3、5、10项的和，观察逼近效果：

import numpy as np import matplotlib.pyplot as plt # 将符号表达式转换为数值函数 f_series = [lambdify(x, fourier_series(N), 'numpy') for N in [1, 3, 5, 10]] # 生成采样点 x_vals = np.linspace(0, 4*np.pi, 1000) # 绘制结果 plt.figure(figsize=(10, 6)) for i, f in enumerate(f_series): plt.plot(x_vals, f(x_vals), label=f'前{2*i+1}项和') # 绘制理想方波 ideal = np.where((x_vals % (2*np.pi)) < np.pi, V, -V) plt.plot(x_vals, ideal, 'k--', label='理想方波') plt.legend() plt.xlabel('x') plt.ylabel('幅值') plt.title('方波的傅里叶级数逼近') plt.grid(True) plt.show()

5. 工程应用与扩展

5.1 PWM波形分析

在电力电子中，PWM（脉宽调制）波形可以看作是方波的变种。我们可以修改之前的代码来分析占空比可变的PWM波：

duty = symbols('d', positive=True) # 占空比(0<d<1) pwm_wave = Piecewise( (V, (x >= 0) & (x < 2*pi*d)), (-V, (x >= 2*pi*d) & (x < 2*pi)), (pwm_wave.subs(x, x - 2*pi), True) ) # 计算PWM波的傅里叶系数 b_n_pwm = (1/pi) * integrate(pwm_wave * sin(n*x), (x, 0, 2*pi)) b_n_pwm = simplify(b_n_pwm.subs(cos(2*pi*d*n), cos(2*d*n*pi)))

5.2 谐波失真分析

通过傅里叶级数，我们可以计算总谐波失真(THD)：

def calculate_thd(N): fundamental = b_n_final.subs(n, 1) harmonics = sum([(b_n_final.subs(n, 2*k-1)/fundamental)**2 for k in range(2, N+1)]) return sqrt(harmonics)

5.3 性能优化技巧

当处理高阶谐波时，符号计算可能变慢。我们可以：

1. 预先计算并存储系数
2. 使用数值积分替代符号积分
3. 并行计算不同谐波分量

from sympy.utilities.lambdify import lambdify from joblib import Parallel, delayed def compute_harmonic(k): n_val = 2*k - 1 coeff = float(b_n_final.subs(n, n_val)) return coeff * np.sin(n_val * x_vals) # 并行计算 results = Parallel(n_jobs=4)(delayed(compute_harmonic)(k) for k in range(1, 21)) waveform = sum(results)

6. 完整代码实现

以下是整合后的完整代码，包含所有功能：

import sympy as sp import numpy as np import matplotlib.pyplot as plt from sympy.utilities.lambdify import lambdify # 符号定义 x, n = sp.symbols('x n', real=True) V = sp.symbols('V', positive=True) # 定义方波 square_wave = sp.Piecewise( (V, (x >= 0) & (x < sp.pi)), (-V, (x >= sp.pi) & (x < 2*sp.pi)), (square_wave.subs(x, x - 2*sp.pi), True) ) # 计算傅里叶系数 b_n = (1/sp.pi) * sp.integrate(square_wave * sp.sin(n*x), (x, 0, 2*sp.pi)) b_n = sp.simplify(b_n.subs(sp.cos(sp.pi*n), (-1)**n)) b_n_final = sp.Piecewise( (0, sp.Eq(n%2, 0)), (4*V/(sp.pi*n), True) ) # 傅里叶级数函数 def fourier_series(N): series = 0 for k in range(1, N+1): series += b_n_final.subs(n, 2*k-1) * sp.sin((2*k-1)*x) return series # 可视化 x_vals = np.linspace(0, 4*np.pi, 1000) V_val = 1.0 # 设幅值为1 plt.figure(figsize=(12, 8)) for N in [1, 3, 5, 10, 20]: fs = fourier_series(N) f_numeric = lambdify(x, fs.subs(V, V_val), 'numpy') plt.plot(x_vals, f_numeric(x_vals), label=f'N={N}') # 理想方波 ideal = np.where((x_vals % (2*np.pi)) < np.pi, V_val, -V_val) plt.plot(x_vals, ideal, 'k--', linewidth=1.5, label='理想方波') plt.title('方波的傅里叶级数逼近', fontsize=14) plt.xlabel('x', fontsize=12) plt.ylabel('幅值', fontsize=12) plt.legend(fontsize=10) plt.grid(True) plt.tight_layout() plt.show()

注意：实际应用中，可以根据需要调整V的值和观察的周期数。代码中的N表示包含的谐波数量，N越大逼近效果越好，但计算量也会增加。

通过这个完整的示例，我们不仅实现了方波傅里叶级数的符号推导，还创建了一个可交互的工具，可以直观地观察不同谐波数量下的逼近效果。这种方法比单纯记忆公式更有助于深入理解傅里叶分析的原理。