news 2026/11/3 17:24:35

用Python复刻傅里叶级数动画:旋转圆叠加如何画出方波

作者头像

张小明

前端开发工程师

1.2k 24
文章封面图
用Python复刻傅里叶级数动画:旋转圆叠加如何画出方波

如果你刷过数学科普类的视频,大概率逃不开那个经典画面:最里面一个小圆在转,圆上又套着第二个圆,再套第三个……好几个圆像机械表盘一样层层咬合,最外层笔尖扫出来的轨迹,从缓慢起伏变得越来越像一条方波。我第一次看到这个动图时愣了很久,脑子里只有一个念头:这到底是怎么画出来的?那会儿刚学完傅里叶级数,公式会背,但根本不知道它"画"的到底是什么。后来为了给学生做演示,我花了一个下午用Python把它从零复刻了一遍,才真正搞懂那些旋转圆和方波之间的关系。这篇文章就把整个复刻过程完整写下来——数学上怎么拆、代码怎么搭、参数怎么调、有哪些文档里不会写的坑,以及最后方波在电路和电机里到底有什么实际关系。无论你是要做教学演示、写可视化程序,还是单纯想把那个经典动画亲手做出来,这篇都可以直接参考。

1. 先看懂画面:那些圆到底在"叠"什么

1.1 画面里其实有两条信息流

这个动画通常分成左右两块。左边是一串连环圈,圆心一个套一个,从原点出发,每个圆都有自己的半径和转速,最外层有一个"笔尖",也就是最后一根向量的端点。右边是一条正在被"画"出来的曲线,笔尖每走一步,曲线就往右延长一点,最终形成一条完整的方波。

我第一次看的时候只注意到那堆圆在转,完全没意识到右边那条方波和圆是一一对应的。后来才明白,左边每一个圆,其实都代表一个正弦分量,也就是方波傅里叶级数里的一项。圆的半径代表这个分量的幅度,圆的转速代表这个分量的频率,整条圆链就是把这些分量"加"在一起的机械结构。换句话说,动画不是"画了一条方波",而是"用一堆圆把方波拼了出来",拼的过程被可视化了而已。

1.2 圆周运动怎么变成正弦波:一句话原理

一条正弦波可以这样理解:让一个点以恒定角速度绕圆圈转,同时把这个点的位置投影到竖直方向,投影点做的往复运动就是正弦波。反过来说,一个正弦波,就是"某个匀速旋转的向量在某一根轴上的影子"。

所以当你把多个不同半径、不同转速的圆按首尾相接的方式叠在一起,每个圆末端的位置就是所有旋转向量的矢量和。矢量和在竖直方向的分量,等于所有正弦分量加起来的值。这就是傅里叶合成(Fourier synthesis)的几何含义:多个旋转向量首尾相接,等价于多个正弦波相加。动画右边那条方波,其实就是整条圆链最后那个笔尖在竖直方向上的坐标随时间变化的轨迹。搞懂这一层之后,剩下的问题就只有一个:方波到底需要哪些"圆",也就是需要哪些分量?

2. 方波的傅里叶级数:公式拆开看只有三件事

2.1 公式先摆出来

一个幅值为±1、周期为2π的理想方波,它的傅里叶级数是:

f(t) = (4/π) × [ sin(t) + (1/3)·sin(3t) + (1/5)·sin(5t) + (1/7)·sin(7t) + ... ]

对应到动画里就是:第一个圆的半径是 4/π,以角速度1转;第二个圆的半径是 4/(3π),以角速度3转;第三个圆的半径是 4/(5π),以角速度5转……依此类推,永远只出现奇数次谐波。

这个公式我第一次看的时候觉得完全是天降公式,但它的每一步其实都有明确逻辑,而且每个数字都有几何含义。你以为要背一堆东西,其实整个公式就围绕三件事:频率取哪些、幅度怎么算、为什么没有偶次项。下面一个个拆。

2.2 为什么只有奇数倍频率

方波有个非常明显的对称性:把波形平移半个周期,得到的新波形等于原波形上下翻转,也就是 f(t+π) = -f(t)。这在信号处理里叫半波反对称,也叫奇谐对称。

现在假设某个偶次谐波存在,它在平移半周期后会变成什么样?sin(n(t+π)) = sin(nt + nπ)。当 n 是偶数时,nπ 是 2π 的整数倍,所以 sin(nt + nπ) = sin(nt),波形平移后完全没变化;当 n 是奇数时,sin(nt + nπ) = -sin(nt),波形平移后正好反相。方波平移半周期要求所有分量都反相,所以只有奇次谐波能满足这个条件,偶次谐波的系数只能全部归零。这个结论完全不需要算积分,就是纯看对称性,逻辑非常干净。也因此你在动画里永远看不到2倍频、4倍频的圆,它们压根没有存在的资格。

2.3 系数4/(nπ) 是怎么来的

再往深一层,每个谐波的幅度系数 4/(nπ) 来自标准的傅里叶系数公式。对周期函数 f(t),正弦项的系数是:

b_n = (1/π) × ∫从 -π 到 π f(t)·sin(nt) dt

把方波分段代入,在 0 到 π 上 f(t)=1,在 -π 到 0 上 f(t)=-1,分段积分后就会得到:n 为偶数时 b_n=0,n 为奇数时 b_n = 4/(nπ)。整个过程就是"分段积分+代入上下限",二十分钟能自己导一遍。导完你会发现一个特别重要的规律:谐波次数越高,幅度按 1/n 衰减。这就是为什么动画里第一个圆最大、后面的圆越来越小,也是为什么你只需要四五个圆就能看出方波的大致形状——高阶谐波只是给方波"修边"而已。

3. 复平面视角:把"圆链"变成程序能执行的逻辑

3.1 旋转向量和复数指数是一回事

要在代码里画旋转圆,最自然的方式是借用复数。一个长度为 r、绕原点以角速度 nω 旋转的向量,在时刻 t 可以写成复数形式:

z = r · e^(i·n·ω·t)

它的实部是 r·cos(nωt),虚部是 r·sin(nωt)。你把虚部当作竖直方向,每个旋转向量在竖直方向上的投影,正好就是傅里叶级数里的 sine 项。整条圆链末端的位置,就是所有向量做复数加法。

所以"方波傅里叶级数动画"本质上就是在画一个复平面,播放从 t=0 到 t=2π 的整个过程,每个时刻把下面这个和式算出来:

S(t) = Σ_{k=0}^{N-1} (4/((2k+1)π)) · e^(i·(2k+1)·t)

S(t) 的虚部就是当前时刻波形的高度,实部则决定了笔尖在圆链里"挂"在什么水平位置。

3.2 圆链的构造规则,其实只有四条

把上面这层逻辑翻译成"画圆"的规则,非常清楚:

  1. 从原点出发,第 k 个圆的圆心放在前 k 个向量的和 S_k(t) 的位置上;
  2. 第 k 个圆的半径 r_k = 4/((2k+1)π),只与谐波序号有关,与时间无关;
  3. 第 k 个圆的旋转角速度是 (2k+1)·t,所以每一帧都要重新计算角度;
  4. 最后一根向量末端 S_{N-1}(t) 的位置就是"笔尖",取它的虚部就是波形纵坐标。

规则比想象中简单。难点只在于:同一时刻,所有圆的圆心位置都是动态的——第二个圆的圆心是第一个向量当前的末端,第三个圆的圆心是前两个向量当前的矢量和……所以画的时候不能一帧一帧硬编码坐标,得让程序从第一个圆开始逐层往下"接力"计算。这也是写代码时最容易犯糊涂的地方,我在第4节里会专门点出那个顺序坑。

4. Python落地:一套能直接跑起来的动画代码

4.1 选型:为什么用 matplotlib 而不是别家

做这种动画,工具选择其实很开放。可以用 JavaScript 加 Canvas 做网页版,也可以用 Processing,但我在教学场景里还是推荐 Python 搭配 matplotlib 的 FuncAnimation:环境几乎人人都有,代码可读性好,画出的矢量图能直接导出 GIF 或视频,要嵌入 PPT 或者公众号都方便。下面这套代码我在 Jupyter 和纯脚本模式下都跑过,可以直接照抄。

需要提前装好 numpy 和 matplotlib;如果想把动画保存成 GIF,还需要 pillow:

pip install numpy matplotlib pillow

4.2 核心函数:一帧状态全在这里

整个程序的灵魂是一个计算函数,输入一个时刻 t 和要用的谐波个数,输出这一时刻所有旋转向量和所有圆心位置:

import numpy as np def epicycle(t, n_harmonics): """给定时刻 t,返回旋转向量列表和圆心位置列表""" vectors = [] # 每一根旋转向量(用复数表示) centers = [0j] # 圆心位置,第一个圆在原点 pos = 0j # 当前累计位置 for k in range(n_harmonics): n = 2 * k + 1 # 只用奇次谐波 amp = 4.0 / (n * np.pi) # 半径 = 振幅 vec = amp * np.exp(1j * n * t) vectors.append(vec) pos += vec # 累加:作为下一个圆的圆心 centers.append(pos) return vectors, centers

这个函数里我最开始写错过的地方是centers.append(pos)和pos += vec的顺序。必须先完成累加再把新位置追加为圆心,顺序一颠倒,整条链就错位一个环节,画出来第二圈圆就跟第一圈脱节了。这类小坑不跑一次动画根本发现不了。

4.3 波形预计算与画布搭建

动画右边要画"已经走过的波形",所以需要一条足够密集的完整曲线,动画播放时按进度截取显示。这里用TRACE_N=1500个采样点,画出一整个周期的波形:

# 全局参数 N_HARMONICS = 12 # 圆的个数 N_FRAMES = 360 # 每周期多少帧 TRACE_N = 1500 # 波形曲线采样密度 # 预计算波形:每个采样点都做一次傅里叶合成 t_dense = np.linspace(0, 2*np.pi, TRACE_N) wave_y = [] for tt in t_dense: _, centers = epicycle(tt, N_HARMONICS) wave_y.append(centers[-1].imag) wave_y = np.array(wave_y) # 动画帧对应的时间点 t_frames = np.linspace(0, 2*np.pi, N_FRAMES)

这里有个新手很容易踩的坑:不要一上来就尝试对整条时间轴做向量化再取虚部。虽然 numpy 能算得更快,但逐点计算更符合"每个时刻独立做一次傅里叶合成"的逻辑,而且这个计算量在十几个圆、一千多个采样点的规模下毫秒级就出结果,完全没有性能压力。先把逻辑写清楚,再谈优化。

4.4 动画更新函数与完整脚本

动画核心是animate(frame),每帧做两件事:左边清空重画圆链,右边按进度画波形。完整代码如下,参数都在文件开头,改N_HARMONICS就能直观看到"圆越多方波越方"的效果:

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation N_HARMONICS = 12 N_FRAMES = 360 TRACE_N = 1500 def epicycle(t, n_harmonics): vectors = [] centers = [0j] pos = 0j for k in range(n_harmonics): n = 2 * k + 1 amp = 4.0 / (n * np.pi) vec = amp * np.exp(1j * n * t) vectors.append(vec) pos += vec centers.append(pos) return vectors, centers t_dense = np.linspace(0, 2*np.pi, TRACE_N) wave_y = np.array([epicycle(tt, N_HARMONICS)[1][-1].imag for tt in t_dense]) t_frames = np.linspace(0, 2*np.pi, N_FRAMES) fig, (ax_left, ax_right) = plt.subplots(1, 2, figsize=(13, 5.5)) plt.subplots_adjust(left=0.08, right=0.95, top=0.88, bottom=0.12) def animate(frame): t = t_frames[frame] vectors, centers = epicycle(t, N_HARMONICS) # ---- 左图:旋转圆链 ---- ax_left.cla() ax_left.set_aspect('equal') # 关键:不设的话圆会变成椭圆 ax_left.set_xlim(-3.4, 3.4) ax_left.set_ylim(-3.2, 3.2) ax_left.set_title(f'{N_HARMONICS} 个圆 / t={t:.2f}') for k in range(N_HARMONICS): r = 4.0 / ((2*k+1) * np.pi) ax_left.add_patch(plt.Circle( (centers[k].real, centers[k].imag), r, fill=False, color='#7d9ad6', lw=1)) v = vectors[k] ax_left.plot([centers[k].real, centers[k+1].real], [centers[k].imag, centers[k+1].imag], color='#f4a261', lw=1.4) ax_left.scatter(0, 0, color='black', s=12) ax_left.scatter(centers[-1].real, centers[-1].imag, color='#e63946', s=28) # ---- 右图:描波形 ---- ax_right.cla() ax_right.set_xlim(0, 2*np.pi) ax_right.set_ylim(-1.5, 1.5) ax_right.set_title('方波逐步成型') ax_right.set_xlabel('t') idx = min(int(frame / N_FRAMES * TRACE_N) + 1, TRACE_N - 1) ax_right.plot(t_dense[:idx], wave_y[:idx], color='#2a9d8f', lw=2) ax_right.scatter(t, centers[-1].imag, color='#e63946', s=28) ax_right.axvline(t, color='#555', ls='--', lw=0.8) anim = FuncAnimation(fig, animate, frames=N_FRAMES, interval=40, blit=False) # 在 Jupyter 里直接 plt.show();保存 GIF 则需要 pillow: # anim.save('fourier_square_wave.gif', writer='pillow', fps=25) plt.show()

这段代码的套路就是每帧cla()清空重画,对 12 个圆来说性能完全够用,逻辑也最容易理解。如果你要跑到 50 个圆以上,建议改成"预先创建图形对象 + set_data 更新"的方式,我在第5节会专门讲。另外在 Jupyter 里跑动画时,记得先把%matplotlib inline换成%matplotlib notebook或%matplotlib tk,不然动画只会显示最后一帧。

4.5 一个重要的坐标解释

刚写完这段代码时,我盯着屏幕疑惑了很久:左边笔尖的横坐标是"所有旋转向量实部之和",右边波形点的横坐标是"时间 t",两个坐标明明不是同一个东西,为什么最终能对上?

答案其实在定义里:左图的横轴是复平面实轴、纵轴是虚轴;右图把时间 t 当作横轴、把笔尖的虚部当作纵轴。两边共用的是"同一时刻 t 下的同一个笔尖位置"——左图强调空间位置,右图强调时间演化。很多第一次做这个动画的人在这个点上绕很久,但只要想明白左右两图共享同一个 t,就通了。如果还不放心,可以把左图笔尖的虚部坐标打出来,和右图当前点的纵坐标对比,两者必然完全相等。

5. 调参指南:让动画既好看又能讲出东西

5.1 谐波数量与吉布斯现象

N_HARMONICS是最值得调的参数。我实测下来:

  • 1个圆:只是一个圆,笔尖画出来是纯正弦;
  • 3个圆:波形开始变"平",边缘还是圆滑的;
  • 7到12个圆:肉眼已经能认出方波,跳变沿明显变陡;
  • 30个圆以上:方波很"方",但如果你把波形放大看,跳变沿附近始终有一圈颤抖的波纹。

这圈波纹就是著名的吉布斯现象:傅里叶级数在间断点附近的过冲不会因为项数增加而消失,极限情况下过冲幅度大约稳定在 9%。很多第一次做这个动画的人以为自己代码写错了,跑来问"为什么圆越多反而有毛刺",其实这是数学本身的性质,不是 bug。这个现象特别适合拿来当教学话题,比单纯看"圆变多波形变方"深刻得多,因为它直接引出了"用有限项逼近一个不连续信号到底能做到多好"的问题。

5.2 布局与视觉细节

有几个让动画"质感"差很多的细节,值得单独说:

  • set_aspect('equal')必须加,否则 matplotlib 会为了让坐标轴填满画布而拉伸,圆会变成椭圆,旋转向量长度看起来也不对,整个数学关系就错乱了;
  • 左图坐标范围要留够。t=0 时所有向量都指向正实轴方向,整条链子会向右伸到接近所有半径之和的位置,12 个圆大概到 2.8 左右,所以我把 x 范围设到 ±3.4,避免动画一开始圆链就跑出画面;
  • 左图的纵轴范围也要配合笔尖的虚部活动范围。笔尖的 y 坐标大约在 ±1 附近,带一点吉布斯过冲,右图的 y 范围设为 ±1.5 就够;
  • 圆的填充建议用透明细线、浅蓝色系,半径臂用橙色,这样圆链的结构一眼能看清;
  • 保存 GIF 时interval和fps要配合好。interval=40表示每帧间隔 40 毫秒,配合fps=25,保存出来的动画节奏和屏显基本一致。

5.3 性能优化:从12个圆到上百个圆

cla()清空重画的好处是代码短、不易出错,坏处是每一帧都要重建所有图形对象。我把N_HARMONICS调到 100 时,动画开始明显掉帧。这时候就要改成"先建好、再更新"的模式:提前创建 N 个Circle对象和 N 条线段,动画函数里只调用set_center、set_data去更新位置,打开blit=True。代码结构稍复杂,但帧率能快一个数量级。

另外,如果只是想要一张静态的方波逼近图,完全不需要动画,直接把N_HARMONICS循环加起来画一条曲线就够了。做这个动画的价值在于"动态"本身——看着圆链转和看着波形画出来,人对谐波叠加的理解完全不是一个层级。

6. 动画之外:方波在现实世界里也有三个舞台

6.1 方波发生电路:信号发生器里的大名

方波不是只在数学动画里出现。实验室里最常见的方波来源,就是用比较器或555定时器搭的方波发生电路:把三角波或正弦波送进比较器,输出端就是陡峭的方波;555在非稳态模式下,靠电容充放电反复翻转输出,也能得到频率可调的方波。

有意思的是,真实电路产生的方波永远不是"理想方波"——上升沿和下降沿总有一定斜率,顶部也可能有轻微振铃。用傅里叶的视角看,边沿越陡,需要的高次谐波就越多。你示波器上测到的那些高频噪声,相当一部分就来自方波自身的谐波成分。所以做信号完整性的人常说,方波是"高频杀手",不是因为它频率高,而是因为它谐波丰富。

6.2 无刷电机的方波调制

另一个我工作里经常打照面的场景,是无刷电机的方波调制。无刷直流电机最简单的驱动方式是六步换相法,每一相在半个电角度周期内通以近似方波的驱动电压,这就是所谓的方波调制。换相处电压突变,电流波形里会携带大量奇次谐波,反映到电机上就是转矩脉动和噪音。

怎么缓解?思路和动画里完全一致:想办法"补"或"滤"掉不需要的高次谐波,或者用正弦波调制代替方波驱动,把谐波成分整体压低。每次我给做电机的同事解释"为什么方波驱动会有噪音"时,直接把傅里叶级数那个旋转圆动画放出来,对方马上就懂了——所有突然跳变,都是无数个圆叠出来的结果。

6.3 为什么工程师也要会"画圆"

最后说点个人体会。做这个动画之前,傅里叶级数对我来说就是课本上"背下来会用"的公式,做工程计算时按部就班套系数。但亲手把这个动画写出来后,我对谐波的理解从一个"抽象记号"变成了一组"可见的向量":每个谐波都有明确的半径、转速和初始相位,叠加的每一帧都看得见摸得着。

这种直观感在工程里非常有用。无论是方波发生电路里选择比较器带宽、无刷电机里评估换相谐波、还是写数字信号处理代码时判断需要保留多少谐波,本质上都在问同一个问题:我需要保留多少个"圆",才能让输出足够好地逼近目标波形?所以有空的话,建议你也花两小时把这个动画写一遍。写完你会发现,后面再遇到和谐波相关的问题,脑子里会多一套"会转的"直觉,而不是一行干巴巴的公式。

版权声明: 本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若内容造成侵权/违法违规/事实不符,请联系邮箱:809451989@qq.com进行投诉反馈,一经查实,立即删除!
网站建设 2026/11/3 17:23:55

解决VSCode中Maven项目Java代码不编译问题

1. 问题现象与背景分析最近在使用VSCode开发Java项目时遇到了一个让人头疼的问题:在一个Maven工程中,src目录下的Java源代码文件无法进行自动编译,也没有代码错误提示(即代码下方不显示红色波浪线)。这意味着即使故意写…

作者头像 李华
网站建设 2026/11/3 17:23:34

Robot Framework接口自动化测试实战与优化

1. Robot Framework 接口自动化测试实战指南作为一名在测试领域摸爬滚打多年的老兵,我见证过太多团队在自动化测试上的挣扎。今天要分享的Robot Framework(后称RF),正是我在多个大型项目中验证过的高效解决方案。不同于市面上那些…

作者头像 李华
网站建设 2026/11/3 17:22:09

专业图床服务与CDN加速技术解析

1. 为什么我们需要专业的图床服务作为一名技术博主,我深知图片加载速度对用户体验的影响有多大。去年我的个人博客就因为图片加载问题流失了近30%的移动端读者。当页面加载时间超过3秒时,53%的用户会选择离开——这是Google的研究数据。传统解决方案通常…

作者头像 李华
网站建设 2026/11/3 17:21:22

深入解析MySQL MVCC机制与实现原理

1. MVCC机制的前世今生第一次接触MVCC这个概念时,我也被各种专业术语绕得头晕。直到某次排查数据一致性问题,亲眼目睹Undo Log如何像时光机一样回溯数据版本,才真正理解这套机制的巧妙设计。MVCC(Multi-Version Concurrency Contr…

作者头像 李华
网站建设 2026/11/3 17:20:50

金属探测器技术解析与选型指南

1. 金属探测器行业现状与选型逻辑 金属探测器作为专业设备,其技术迭代速度远超普通消费者想象。过去五年间,行业平均探测精度提升了47%,而误报率下降了63%。这种进步主要得益于三大技术突破: 1.1 多频同步技术的实际价值 传统单…

作者头像 李华
网站建设 2026/11/3 17:20:13

推理模型与通用模型:CoT、提示语策略与选型指南

简介:《清华大学 DeepSeek从入门到精通》是一份面向AI初学者、产品与技术人员及高校师生的PDF课程材料,围绕DeepSeek公司与DeepSeek-R1开源推理模型展开,帮助读者建立从基础认知到实际应用的完整知识框架。内容系统梳理DeepSeek是什么、能做什…

作者头像 李华