微实践 - 方波信号的傅里叶正弦展开

  法国数学家傅里叶指出任何周期函数都可以表示为不同频率的正弦函数和（或）余弦函数的和，每个正弦或者余弦函数乘以不同的系数。这一看似简单的数学表达被称为傅里叶级数，它构成了数字信号处理的数学基石。

版权声明

本文可以在互联网上自由转载，但必须：注明出处(作者：海洋饼干叔叔)并包含指向本页面的链接。

本文不可以以纸质出版为目的进行改编、摘抄。

  本文节选自《Python编程基础及应用》，高等教育出版社，作者：陈波，刘慧君

  此处我们以数字电路中广泛使用的方波为例，使用不同频率和振幅的正弦波来构造方波。从而帮助读者在高等数学的理论教学之外，体验傅里叶级数的实践价值。

  如图10-4所示的振幅A=1，周期T=2π的标准方波的傅里叶级数展开式为：

  以上述公式为基础，下述程序生成了4个不同频率、不同振幅的正弦信号分量，并将这些分量求和，以“求得”方波信号。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#fourier.py
import numpy as np, matplotlib.pyplot as plt

t = np.linspace(-2*np.pi, 2*np.pi, 1000)
N = 4
r = np.zeros(1000)

plt.xlim(-2*np.pi, 2*np.pi+5)
for n in range(1,N+1):
    f = 4 / np.pi * np.sin((2 * n - 1) * t) / (2 * n - 1)
    r += f
    if n<= 4:
        s = "${4sin(%dt)}/{%d\\pi}$" % (2*n-1,2*n-1)
        plt.plot(t, f+2.0*(5-n),label=s)

plt.plot(t, r, label='$r=\sum$')
plt.legend()
plt.title("Fourier Series of a Square Wave")
plt.show()

  上述程序的执行结果如图10-5所示。

🚩第4行：生成从-2π到+2π，包含1000个元素的等差数列t。数组t在本例中对应时间轴（横轴）。

🚩第5行：N表示拟参与叠加生成“方波”的正弦分量数。

🚩第6行：生成形状为(1000,)的浮点数数组r，用于保存结果“方波”。

🚩第8行：设置matplotlib图表的横轴范围。在本例中，时间轴t的取值范围为-2π ~ +2π，请注意本行代码故意给右限值加5，是为了给图10-5右上角的图示框预留空间。

🚩第9 ~ 11行：通过循环依次计算各正弦分量，并累加给“方波”r。

🚩第10行：按公式计算方波的正弦分量f。根据10.1.4节中的讨论，本行中的np.sin()函数，以及乘、除等四则运算，均为ufunc函数，其作用范围为数组的全部元素。

🚩第11行：将正弦分量f累加给“方波”r。

🚩第12 ~ 14行：绘制前4个正弦分量曲线。

🚩第13行：s对应分量曲线的标签。请注意，该标签字符串以美元符号开始，美元符号结束，这表示该标签字符串为LaTeX格式公式。LaTeX公式是在学术文献写作过程中广泛采用的一种数学公式标记方法，详情可网上查询。由于LaTeX公式采用{}来包裹公式符号，正好与4.1.6节所述的字符串格式化中所用占位符形式冲突。为了处理这种冲突，我们不得不采用%格式的字符串格式化方法，语法细节请见4.1.6节中的扩展阅读。

🚩第14行：plt.plot(t, f+2.0*(5-n),…)函数用于绘制折线图，当折线图的构成点足够密集时，看起来就是一条光滑曲线。t对应所有离散点的横轴坐标，f+2.0*(5-n)则对应所有离散点的纵轴坐标，关键字参数label为曲线的标签。为了避免各曲线在图表中重叠显示，这里给正弦分量f加上了2.0*(5-n)，有意将相应曲线在纵向上抬高。请读者注意，这里的f上形状为(1000,)的一维数组，而2.0*(5-n)的结果为浮点数，两者之间的加法操作事实上是通过ufunc函数进行的，该加法函数将f中的每个元素都与浮点数相加，生成一个新的一维数组。

🚩第16行：以标签“r=∑”显示结果“方波”r。

  如图10-5中的结果方波“r=∑”看起来并不标准，因为它只包含了4个不同频率和振幅的正弦分量。按傅里叶级数，参与叠加的正弦分量的数量增多，结果曲线越接近标准方波。图10-6展示了包含20个和1000个正弦分量的方波信号。在程序fourier.py中，修改N值即可调整方波所包含的正弦分量个数。