别再死记硬背公式了！用Python的NumPy和Matplotlib亲手‘画’出傅里叶级数（附完整代码）

用Python动态可视化傅里叶级数：从数学公式到交互式图形

记得第一次接触傅里叶级数时，那些复杂的公式让我头晕目眩——直到我亲手用代码将它"画"出来。本文将带你用Python的NumPy和Matplotlib，通过编写可交互的代码，直观理解这个强大的数学工具如何将任意周期函数分解为简单的正弦余弦波组合。

1. 准备工作：搭建Python分析环境

在开始之前，我们需要配置好Python环境并安装必要的库。推荐使用Jupyter Notebook或Google Colab这类交互式环境，可以实时看到代码执行结果。

核心工具包安装：

pip install numpy matplotlib ipywidgets

表：本教程使用的主要Python库及作用

import numpy as np import matplotlib.pyplot as plt from ipywidgets import interact, IntSlider %matplotlib inline

提示：如果你在本地运行遇到显示问题，可以尝试在代码开头添加%matplotlib notebook以获得更好的交互体验。

2. 从方波开始：理解傅里叶级数的基本原理

让我们从一个经典的例子入手——方波的傅里叶级数展开。方波在信号处理中非常常见，其数学表达式为：

def square_wave(t, T=2*np.pi): """生成周期为T的方波""" return np.where(np.sin(2*np.pi*t/T) > 0, 1, -1)

方波的傅里叶级数展开式为： $$ f(t) = \frac{4}{\pi}\left[\sin(\omega t) + \frac{1}{3}\sin(3\omega t) + \frac{1}{5}\sin(5\omega t) + \cdots\right] $$

其中ω=2π/T是基频。让我们用Python实现这个级数的前N项和：

def fourier_series_square(t, N, T=2*np.pi): """计算方波傅里叶级数的前N项和""" omega = 2*np.pi/T result = np.zeros_like(t) for n in range(1, N+1, 2): # 只包含奇数次谐波 result += (4/(n*np.pi)) * np.sin(n*omega*t) return result

3. 动态可视化：观察级数如何逼近原函数

现在我们可以创建一个交互式可视化，观察随着项数N的增加，傅里叶级数如何逐步逼近方波：

def plot_fourier_approximation(N=5): t = np.linspace(-2*np.pi, 2*np.pi, 1000) original = square_wave(t) approximation = fourier_series_square(t, N) plt.figure(figsize=(10, 6)) plt.plot(t, original, 'r-', linewidth=2, label='原始方波') plt.plot(t, approximation, 'b-', label=f'傅里叶级数(N={N})') plt.title('方波的傅里叶级数逼近') plt.xlabel('时间 t') plt.ylabel('幅值') plt.legend() plt.grid(True) plt.show() interact(plot_fourier_approximation, N=IntSlider(min=1, max=101, step=2, value=1))

运行这段代码，你会看到一个滑块控件，拖动它可以实时观察不同项数下的逼近效果。当N=1时，我们只有一个简单的正弦波；随着N增加，越来越多的谐波被加入，波形逐渐接近理想的方波。

4. 深入原理：分解与合成的数学实现

傅里叶级数的核心思想是将复杂函数分解为不同频率的正弦余弦波。让我们看看如何计算任意周期函数的傅里叶系数：

def compute_fourier_coeffs(func, N, T=2*np.pi): """计算周期函数func的前N个傅里叶系数""" omega = 2*np.pi/T a_coeffs = [] b_coeffs = [] # 计算a0 t = np.linspace(0, T, 1000) a0 = 2/T * np.trapz(func(t), t) a_coeffs.append(a0) # 计算an和bn for n in range(1, N+1): an = 2/T * np.trapz(func(t)*np.cos(n*omega*t), t) bn = 2/T * np.trapz(func(t)*np.sin(n*omega*t), t) a_coeffs.append(an) b_coeffs.append(bn) return a_coeffs, b_coeffs

我们可以用这个函数分析任意周期信号。例如，分析一个三角波：

def triangle_wave(t, T=2*np.pi): """生成周期为T的三角波""" return 2*np.arcsin(np.sin(2*np.pi*t/T))/np.pi # 计算前10个傅里叶系数 a_coeffs, b_coeffs = compute_fourier_coeffs(triangle_wave, 10) # 打印结果 print("傅里叶系数a_n:", a_coeffs) print("傅里叶系数b_n:", b_coeffs)

5. 扩展应用：分析真实世界信号

傅里叶分析不仅适用于理想波形，也可以用于分析真实信号。让我们模拟一个包含噪声的信号并分析其频率成分：

def analyze_real_signal(): # 生成含噪声的复合信号 t = np.linspace(0, 2*np.pi, 1000) signal = (np.sin(t) + 0.5*np.sin(3*t) + 0.3*np.sin(7*t) + np.random.normal(0, 0.2, len(t))) # 计算傅里叶系数 N = 20 a_coeffs, b_coeffs = compute_fourier_coeffs(lambda x: signal, N) # 绘制结果 plt.figure(figsize=(12, 8)) plt.subplot(2, 1, 1) plt.plot(t, signal, label='原始信号') plt.title('含噪声的时域信号') plt.xlabel('时间') plt.ylabel('幅值') plt.legend() plt.subplot(2, 1, 2) frequencies = np.arange(N+1) plt.stem(frequencies, a_coeffs, 'r', markerfmt='ro', label='a_n (余弦系数)') plt.stem(frequencies[1:], b_coeffs, 'b', markerfmt='bo', label='b_n (正弦系数)') plt.title('傅里叶系数频谱') plt.xlabel('谐波次数 n') plt.ylabel('系数值') plt.legend() plt.tight_layout() plt.show() analyze_real_signal()

这个例子展示了如何从嘈杂的信号中提取出其主要频率成分。通过分析傅里叶系数，我们可以识别出信号中存在的各个频率分量及其相对强度。

6. 性能优化：使用FFT加速计算

对于长信号序列，直接计算傅里叶系数效率较低。NumPy提供了快速傅里叶变换(FFT)算法来加速这一过程：

def fft_analysis(signal, sample_rate): """使用FFT进行频谱分析""" n = len(signal) fft_result = np.fft.fft(signal)/n freqs = np.fft.fftfreq(n, 1/sample_rate) # 只取正频率部分 half_n = n//2 return freqs[:half_n], np.abs(fft_result[:half_n]) # 生成测试信号 sample_rate = 1000 # 采样率1kHz t = np.linspace(0, 1, sample_rate, endpoint=False) signal = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t) # 进行FFT分析 freqs, amplitudes = fft_analysis(signal, sample_rate) # 绘制频谱图 plt.figure(figsize=(10, 4)) plt.plot(freqs, amplitudes) plt.title('FFT频谱分析') plt.xlabel('频率 (Hz)') plt.ylabel('幅值') plt.xlim(0, 200) plt.grid() plt.show()

这段代码展示了如何使用FFT快速分析信号的频率成分。在实际工程应用中，FFT是处理信号频谱的标准工具。

7. 高级应用：交互式傅里叶合成器

最后，我们创建一个更复杂的交互式工具，允许用户自定义各谐波的权重，实时观察合成效果：

from ipywidgets import FloatSlider, VBox def interactive_synthesizer(): # 创建滑块控件 sliders = [] for n in range(1, 11): slider = FloatSlider(min=0, max=1, step=0.05, value=0, description=f'Harmonic {n}', continuous_update=True) sliders.append(slider) # 合成函数 def synthesize(**kwargs): t = np.linspace(-np.pi, np.pi, 1000) result = np.zeros_like(t) for n, weight in kwargs.items(): harmonic_num = int(n.split('_')[1]) result += weight * np.sin(harmonic_num * t) plt.figure(figsize=(10, 5)) plt.plot(t, result) plt.title('傅里叶合成结果') plt.xlabel('时间') plt.ylabel('幅值') plt.grid() plt.ylim(-5, 5) plt.show() # 创建交互界面 interact(synthesize, **{f'h_{i}': slider for i, slider in enumerate(sliders, 1)}) interactive_synthesizer()

这个工具可以让你直观地体验傅里叶合成的过程，通过调整不同谐波的权重，创造出各种复杂的波形。