☰
用Python替代MATLAB完成信号与系统时域分析
2026/10/14 19:51:25 网站建设 项目流程

学信号与系统这门课,时域分析是第一个绕不开的坎。每逢实验周,最常见的情况是大家齐刷刷打开 MATLAB,结果要么机器上没装授权,要么大作业提交前才发现导出的图片惨不忍睹。其实用 Python 加 matplotlib 就能把这一套实验做得很顺手:波形生成、时移反转、卷积运算、LTI 系统零状态响应,全都能在几十行代码内完成,而且免费、跨平台、方便复盘。这篇文章不是课程讲义,而是我从完成一次“信号的时域分析”实验里整理出的完整实操记录,覆盖环境搭建、典型信号生成、时域运算、卷积与系统响应,以及大量课本上不会写的细节。无论你是正在做实验抓耳挠腮的同学,还是想快速回顾时域分析工程实现的人,都可以照着跑一遍,再换成自己的参数去试。

1. 时域分析到底要做什么

1.1 先理清时域分析的核心脉络

信号与系统里的“时域”,字面意思就是以时间为自变量去看信号和系统。输入信号用 x(t) 表示,输出用 y(t) 表示,系统的作用被抽象成一种对时间波形的“加工”。时域分析的本质,就是直接研究这种加工过程:信号长什么样、做了哪些运算、系统冲激响应如何、零状态响应怎么算。听起来很理论,但一旦落到图和代码里就特别直观。

在这个实验里,我习惯把任务拆成四层。第一层是典型信号生成,比如单位阶跃、单位冲激、正弦、指数衰减、矩形脉冲,这些是最基础的构建块。第二层是信号的基本运算,时移、反转、尺度变换、加减乘,目的是理解数学表达式和波形之间的一一对应关系。第三层是卷积,这是时域分析的重中之重,两个信号相卷积的结果,经常就是某个系统对输入的响应。第四层是线性时不变系统分析,给定系统冲激响应 h(t),求输入 x(t) 作用下的零状态响应 y(t),卷积的计算价值在这里体现得最明显。

很多同学做实验时容易陷入“只要画出图就算完成”的状态。其实画图只是验证手段,真正要掌握的是怎么把连续时间数学公式变成计算机里可处理的数组。这里存在一个关键的思维切换:计算机里没有真正的连续信号,只有密集采样后的离散序列。所以做时域分析实验,本质上是在用离散数值方法近似连续时间问题,理解这一点,后面很多参数设置就顺理成章了。

1.2 实验任务拆解与工具链选择

具体到实验报告,常见要求无非这几项:

  • 画出常见连续时间信号和离散时间信号的波形。
  • 验证信号时移、反转、尺度变换对波形的影响。
  • 计算信号的能量与平均功率。
  • 用卷积求两个信号的卷积结果,画出卷积波形。
  • 给定系统冲激响应或微分方程,求输入信号作用下的零状态响应。

每一项用 Python 都有对应的库。NumPy 负责数组运算和信号采样,SciPy 提供卷积、滤波、线性系统仿真等高级工具,matplotlib 负责出图。用过 MATLAB 的人会发现映射关系很直接:MATLAB 的plot对应plt.plot,stem对应plt.stem,conv对应np.convolve,lsim对应scipy.signal.lsim。基本不需要额外学什么复杂的语法,上手阻力很小。

工具选型上我坚持用 Python,主要有三个原因。一是免费,学生自己电脑上没有授权限制,实验室批量安装也容易;二是过程透明,每一步中间变量都是普通数组,你可以随时print出来检查,不像黑盒工具箱那样难以调试;三是复用性好,同样的代码稍微改改就能继续做后续的频谱分析、离散系统实验,不需要换语言重新学习一套绘图接口。当然,MATLAB 在专业工具箱和文档完整性上仍然有优势,但针对本科阶段的时域分析实验,Python 已经绰绰有余。

2. 环境准备与典型信号生成

2.1 环境搭建与连续信号仿真思路

动手写代码前,先把环境弄干净。我用的是 Python 3.9 以上版本,依赖就三个:numpy、matplotlib、scipy。安装命令很常规:

pip install numpy matplotlib scipy

如果你用的是 Anaconda 这类集成环境,这三个包一般已经预装。版本方面不用追新,numpy 1.24 以上、matplotlib 3.7 以上、scipy 1.10 以上,足够覆盖实验需求。

连续信号仿真的思路,我在前面已经点破了:用密集的离散时间点近似连续时间。比如想模拟从 0 秒到 10 秒的时间范围,可以这样定义时间轴:

import numpy as np dt = 0.01 # 采样间隔,越小波形越平滑 t = np.arange(0, 10 + dt, dt)

这里dt是仿真采样间隔,可以把它看成相邻两个采样点之间的时间步长。dt越小,曲线越光滑,但计算量也越大。对于普通实验波形,dt = 0.01已经是很好的精度,画出来人眼看不出折线感。如果要模拟更高频的信号,就得相应减小dt,保证一个周期里至少有 10 个以上采样点,否则波形会明显失真。

有一个很容易被忽略的细节:用np.arange(start, stop, dt)时,stop通常不会刚好落在你想要的终点上。我用10 + dt作为终点,就是为了让最后一个点尽量接近 10。如果直接写成np.arange(0, 10, 0.01),最后一个点可能是 9.99,虽然差别极小,但计算信号周期或做积分时,端点差异会带来少量误差。

2.2 典型信号的 Python 实现

各类典型信号是时域分析的基石。我在这里给出最常用的一组写法,你可以直接写成一个函数库,后续所有实验都复用。

单位阶跃信号最容易写,也有很多种表达:

u = (t >= 0).astype(float) # t>=0 为 True,赋值成 1.0 # 或者用 numpy 自带函数 u2 = np.heaviside(t, 0.5) # 在 t=0 处取 0.5

两种写法在多数场景下等价,只有 t=0 这点不同。严格讲,连续单位阶跃信号在原点的定义有多种约定,常见教材取 1,也有取 0.5 的。做数值仿真时,(t >= 0)这种写法让 t=0 处为 1,适合大多数时域响应计算;np.heaviside(t, 0.5)则更符合某些定义。到底用哪个,取决于你的教材约定,但不要两种混用,否则报告里会出现莫名其妙的 0.5 偏差。

矩形脉冲可以基于阶跃信号构造,也可以用布尔索引:

rect = ((t >= 0) & (t < 2)).astype(float) # 0~2 秒内为 1

正弦信号和指数衰减信号则更直接:

A, f, phi = 1.0, 1.0, 0.0 x_sin = A * np.sin(2 * np.pi * f * t + phi) alpha = 1.0 x_exp = np.exp(-alpha * t) * (t >= 0)

冲激信号是时域实验里最常出问题的一个。连续时间单位冲激 δ(t) 是一个理想化信号,计算机里无法直接生成,常规做法是用一个宽度为dt、高度为1/dt的窄矩形脉冲去近似它。这样做的目的是保证脉冲面积为 1:1/dt * dt = 1。代码长这样:

delta = np.zeros_like(t) delta[np.abs(t) <= dt / 2] = 1.0 / dt

画出来是一根近似“针”的矩形。如果你把dt取得很小,高度会非常大,这容易让人误以为是程序 bug。其实这是正常现象,连续冲激的“强度”体现在面积而不是某一点的高度,近似脉冲的高度由采样间隔决定,dt越小越接近理想冲激。

抽样信号sinc在时域分析里也常用,用 numpy 写特别方便:

T = 1.0 s = np.sinc(t / T)

注意np.sinc(x)定义为sin(pi*x)/(pi*x),和通信领域的Sa(x)=sin(x)/x差了一个 π 因子。用之前务必看公式对应的定义,否则画出来的曲线位置会跟教材对不上。

2.3 画图细节:连续曲线与离散杆图

信号生成后,接下来就是画图。matplotlib 的plot适合连续时间信号,画出来是光滑折线;离散时间信号则要画成离散序列,用stem更专业。

以一组典型信号对比图为例:

import matplotlib.pyplot as plt plt.rcParams['font.sans-serif'] = ['SimHei'] # 防止中文乱码,按系统调整 plt.rcParams['axes.unicode_minus'] = False # 负号正常显示 fig, axs = plt.subplots(2, 3, figsize=(12, 6)) axs[0, 0].plot(t, u, linewidth=1.2) axs[0, 0].set_title('单位阶跃信号') axs[0, 1].plot(t, rect, linewidth=1.2) axs[0, 1].set_title('矩形脉冲信号') axs[0, 2].plot(t, x_sin, linewidth=1.2) axs[0, 2].set_title('正弦信号') axs[1, 0].plot(t, x_exp, linewidth=1.2) axs[1, 0].set_title('指数衰减信号') axs[1, 1].plot(t, delta, linewidth=1.2) axs[1, 1].set_title('冲激信号近似') # 离散信号示例 n = np.arange(-5, 6) xn = np.cos(0.5 * n) axs[1, 2].stem(n, xn, linefmt='-', markerfmt='o', basefmt=' ') axs[1, 2].set_title('离散时间余弦信号') for ax in axs.flat: ax.grid(True, linestyle='--', alpha=0.6) ax.set_xlabel('t/s') plt.tight_layout() plt.show()

画离散信号时,plt.stem的三个参数值得留意:linefmt控制竖线样式,markerfmt控制圆点样式,basefmt控制基线。基线的默认样式是一条实线,但隔很远看就很碍眼,所以用basefmt=' '去掉。新版 matplotlib 的stem已经不需要use_line_collection=True这个历史参数,加了反而可能收到警告,直接删掉即可。

绘图规范看起来是小事,但实验报告里的图是否专业,很多时候就差在这些细节。坐标轴标签、网格线、图例、标题一个都不能少,否则老师拿到图根本不知道你画的是什么参数下的信号。

3. 信号的基本时域运算实现与可视化

3.1 时移、反转与尺度变换的向量化实现

时域运算里最基础的是时移、反转和尺度变换。用 Python 实现时,最大的坑在于“数组索引对应时间点”,不能像数学公式那样随意替换自变量。

先说时移。若原信号是 x(t),时移后的信号是 y(t)=x(t-t₀),向右平移时 t₀ 为正。核心思路是:在新时间轴 t 上的每个点,取原信号在 t-t₀ 处的值。如果原信号在同一个离散时间轴 t 上已经定义好,可以直接用np.interp插值:

t0 = 1.0 y_shift = np.interp(t - t0, t, x, left=0, right=0)

np.interp的四个参数依次是查询点、原时间轴、原信号值、越界时的填充值。这里把越界部分填成 0,避免时移出边界时出现意想不到的“卷绕”。如果不写left和right,numpy 会默认用最边上的值填充,导致波形两端出现不该有的平台,这是非常隐蔽的错误。

反转信号的实现同理:

y_reverse = np.interp(-t, t, x, left=0, right=0)

尺度变换 y(t)=x(a·t),关键是理解时间轴压缩和拉伸的方向。a>1 时波形被压缩,a<1 时波形被拉伸。代码还是插值:

a = 2.0 y_scale = np.interp(a * t, t, x, left=0, right=0)

在实际实验报告里,最好把原始信号和运算后的信号画在同一张图上,用不同颜色区分,再在标题里注明参数。比如“若 a=2,波形在时间轴上压缩一半”,这类结论直接看波形就能记住。

3.2 加减乘运算与时间轴对齐

信号相加、相乘看起来简单,实际中却容易因为时间轴没对齐而出错。如果两个信号共用同一个t数组,那直接做数组加减乘就可以:

x1 = np.sin(2 * np.pi * 1.0 * t) x2 = (t >= 0).astype(float) x_add = x1 + x2 x_prod = x1 * x2

但如果你从两个不同的实验里拿到了长度不同的信号,想加在一起,就必须先把它们插值到公共时间轴上。一个稳妥做法是以较密的那个时间轴为基准,用np.interp把另一个信号重采样过去。千万不要直接x1[:len(x2)] + x2,那样会悄悄丢弃采样点,而且时间戳完全错位。

乘法的典型应用是“信号调制”。比如用一个高频余弦信号去乘矩形脉冲,就得到了一个“门控”后的振荡波形:

x_mod = rect * np.cos(2 * np.pi * 10.0 * t)

这类运算在信号与系统后续的频谱搬移章节里会反复出现。时域实验里先看明白波形,等到频域实验里再理解频移规律,会顺很多。

3.3 信号能量与功率的数值计算

时域分析里除了波形,还经常要算能量和平均功率。连续时间信号的能量定义是平方信号在全时域上的积分,数值上用梯形法或简单矩形法都能实现:

E = np.trapezoid(np.abs(x_sin) ** 2, t)

注意 numpy 的接口有过一次调整:老版本是np.trapz,新版本推荐np.trapezoid。如果你用的版本比较老,写成np.trapz也能跑,但会提示 DeprecationWarning。

对于有限能量的非周期信号,平均功率按总能量除以持续时间算。比如前面那个矩形脉冲,如果高度为 1、宽度为 2 秒,那么能量约为 2,平均功率约为 1。对周期信号来说,平均功率要在一个周期上积分再除以周期。以正弦信号 A·sin(2πft) 为例,理论平均功率是 A²/2,有效值是 A/√2。写代码时可以用足够多的整周期样本去验证这个公式,进一步确认数值积分的正确性。

我自己做实验时会习惯先理论手算出结果,再用代码跑一遍做交叉验证。比如正弦幅值 A=2,理论功率就是 2,数字积分结果应该非常接近 2.00。如果差异超过 1%,通常不是计算问题,而是时间轴没有取到整数周期,或者dt太大导致的数值误差。

4. 卷积与系统时域响应分析

4.1 卷积的直观理解与离散实现

卷积是时域分析里最核心的运算。它的数学形式是积分,但如果只看公式,很容易迷失在求和符号里。我更愿意把它理解成“一个信号翻转后,从另一个信号上滑过去,在每个时刻求重叠区域加权面积”。这个比喻虽然不严格,但对理解波形形状非常有帮助。

离散卷积在 numpy 里一行就好:

x = np.array([2.0, 1.0]) h = np.array([1.0, 3.0, 1.0]) y = np.convolve(x, h, mode='full')

输出长度等于len(x) + len(h) - 1,结果是 [2, 7, 5, 1]。为什么是这个长度?因为卷积本质上是在做多项式乘法,两个序列“错位相乘相加”后,横向跨度必然是两个原始跨度之和减 1。

如果只是在课堂上做验证,mode='full'最直观。但如果你想得到跟输入长度一样的输出,可以用mode='same',它会从 full 结果里剪出一段。选取的起点规则有点绕,我的建议是:先用 full 把完整波形看清楚,再考虑要不要 same。

4.2 连续卷积近似:为什么要乘以 dt

连续卷积定义是积分,而计算机只能做求和。把积分中的微元 dt 显式地提出来,这是连续卷积实验里最容易丢分的地方。

例如,宽度为 2 秒、高度为 1 的矩形脉冲与自身卷积,会得到一个三角形脉冲,顶点高度应该是 2。如果你直接写np.convolve(rect, rect),得到的峰值大约是 200(因为两个 rect 各有约 200 个采样点)。只有写成:

y = dt * np.convolve(rect, rect, mode='full')

峰值才会回到 2 附近。这里乘dt的本质,是用离散和近似连续积分时要补上时间步长。忘记这一步,所有卷积结果的纵轴都会整体放大,而且放大倍数还跟采样密度有关,属于最容易暴露“照抄代码却没理解原理”的错误。

除了乘 dt,连续卷积结果的时间轴也要重新计算。如果两个序列都是从 t=0 开始采样,卷积结果的第 0 个点对应 t=0,之后每个点相差 dt,所以时间轴是:

t_conv = np.arange(len(y)) * dt

但如果输入序列的起点不是 0,比如 x 从 t=-2 开始、h 从 t=1 开始,那么卷积输出的起点理论上要相加,而不是默认为 0。对于大多数课程实验,建议大家把信号统一设成从 0 开始,省掉这个坑。

4.3 完整案例:矩形脉冲通过一阶 RC 系统的响应

下面给一个可以直接跑、能写进实验报告的综合案例。假设输入是一个 0 到 2 秒的矩形脉冲,幅度为 1,系统是一阶 RC 低通网络,时间常数 RC=1 秒。这个系统的冲激响应为:

h(t) = (1/RC) * e^(-t/RC) * u(t)

用卷积求零状态响应,再和微分方程数值解做个对比,能同时验证两个方法。代码框架如下:

import numpy as np import matplotlib.pyplot as plt RC = 1.0 dt = 0.01 t = np.arange(0, 12 + dt, dt) # 输入矩形脉冲 x = ((t >= 0) & (t < 2)).astype(float) # 冲激响应 h = (1 / RC) * np.exp(-t / RC) * (t >= 0) # 卷积求响应 y_conv = dt * np.convolve(x, h, mode='full') t_conv = np.arange(len(y_conv)) * dt # 直接用欧拉法解微分方程:RC * y' + y = x y_euler = np.zeros_like(t) for i in range(1, len(t)): dy = (x[i - 1] - y_euler[i - 1]) / RC y_euler[i] = y_euler[i - 1] + dy * dt # 绘图 plt.figure(figsize=(8, 5)) plt.plot(t, x, label='输入 x(t)', linestyle='--', linewidth=1.5) plt.plot(t_conv, y_conv, label='卷积结果', linewidth=1.5) plt.plot(t, y_euler, label='欧拉法数值解', linewidth=1.5) plt.axvline(RC, color='gray', linestyle=':', label='RC 时间常数') plt.xlim(0, 8) plt.ylim(0, 1.2) plt.xlabel('t / s') plt.ylabel('幅度') plt.title('矩形脉冲通过一阶 RC 系统的时域响应') plt.grid(True, linestyle='--', alpha=0.6) plt.legend() plt.tight_layout() plt.show()

这里要特别留意卷积结果只画到 0~8 秒。为什么不全画?因为全np.convolve的结果长度接近 2400 个点,时间轴到了 24 秒,但系统响应在 8 秒之后已经非常接近零,画出来只是给图面右侧加一条长尾巴,反而看不出细节。做实验时你应该主动控制横轴范围,让关注的信息尽量占满画布。

从结果波形上能看到几个典型特征:矩形脉冲的上升沿经过 RC 系统后变成了指数上升,不再陡峭;脉冲结束后输出没有立刻归零,而是按时间常数 RC 指数衰减;峰值大约出现在 2 秒附近,但比输入跳变略有延迟。这些现象都可以在实验报告里定性解释,再和理论冲激响应卷积结果对照。

5. 常见问题与排查技巧实录

5.1 中文乱码与负号显示

新手画图最常遇到的就是中文变方框、坐标轴负号显示出乱码。在代码开头做好全局设置即可:

plt.rcParams['font.sans-serif'] = ['SimHei', 'Noto Sans CJK SC'] plt.rcParams['axes.unicode_minus'] = False

第一行让 matplotlib 支持中文字体,第二行让负号用 ASCII 方式显示。Windows 下一般用SimHei,Linux 或 macOS 可能要用系统中文字体名。如果设置完还是乱码,就换英文标签,实验报告里也能接受,毕竟目标是表达清楚波形含义。

5.2 波形失真与数值坑

数值实验最大的敌人是采样率不足。画正弦信号时,如果dt太大,一个周期里只有 4、5 个点,折线会明显变成“锯齿”,更严重的是视觉上周期可能被看错。我的经验是:一个周期至少采样 20 个点,画面才足够平滑,也基本不会引起误解。对于频率为 f 的信号,时间步长至少要小于 1/(20 f)。

卷积计算慢是另一个高频问题。np.convolve用的是直接法,当两个序列长度都上万时能明显感觉到卡顿。如果实验中需要长序列卷积,换成 FFT 加速版:

from scipy.signal import fftconvolve y = dt * fftconvolve(x, h, mode='full')

结果和直接卷积在数值误差范围内一致,但速度快很多。长序列尤其是做系统长时响应时,这个替换几乎感觉不到差别。

关于连续卷积的幅值和时间轴,我再啰嗦一遍:别忘了乘dt,别忘了重新构造卷积时间轴。这两个问题我见过太多次,几乎成了实验报告的标准扣分点。

5.3 实验报告与出图技巧

实验报告里的图,不仅要“能看”,最好还要“好看且信息完整”。我的常规做法是保存成 PDF 矢量图而不是 PNG 截图,这样在文档里放大也不会糊。保存代码一行就行:

plt.savefig('result.pdf', bbox_inches='tight')

bbox_inches='tight'会自动去掉四周多余空白,让图片紧凑。如果保存 PNG,建议dpi=300。

多信号对比时,用不同线型和透明度区分。实线画结果,虚线画输入,点线画理论参考,图例放在不遮挡曲线位置。必要时用plt.axvline或plt.axhline标记特殊时间点和幅值,比如 RC 时间常数、峰值点。一张图容纳三四条曲线已经接近极限,再多就乱,不如拆成子图。

6. 这些实验做完后的一点体会

6.1 我踩过最值得说的一次坑

第一次跑一阶 RC 系统响应时,我没有乘dt,直接拿矩形脉冲和冲激响应做卷积,结果输出波形幅度大得离谱。当时的第一个念头是“代码写错了”,检查半天却发现np.convolve的调用完全正确。后来才反应过来,连续卷积的离散近似必须要带时间步长,这个问题不是语法错误,而是数学模型到数值模型的映射问题。从那次以后,我每写一步代码都会先在纸上把离散化公式写出来,再对照着填参数。

另一个让我印象深刻的坑是时间轴偏移。离散卷积输出长度等于两序列长度之和减一,这个结论大家都知道,但真正画图时,很多人会草率地沿用输入的时间轴,结果图形右移了一大截。如果你把输入 x(t) 和卷积结果 y(t) 画在同一张图上,横轴必须用各自的正确时间轴,不能偷懒共用一个t。

6.2 这套方案还能怎么扩展

做完时域分析实验之后,后面的很多内容都可以复用这套代码。最简单的扩展是把信号生成函数封装成一个模块,之后做频域分析时继续调用,省得每个实验都从开头矩阵慢慢敲。进阶一点,可以用scipy.signal.lsim替代手写欧拉法,直接仿真任意传递函数的系统响应;也可以用np.correlate或scipy.signal.correlate做相关运算,理解相关和卷积之间的镜像关系。

如果你正在做信号与系统实验,我的建议是不必急着抄完整代码,而是先自己手算一个简单例子的理论值,再用代码验证。这样虽然多花十几分钟,但波形图上的每一条曲线在你眼里都会“透明”起来。后面再遇到奇怪的结果,你也能更快判断是参数问题、采样问题,还是模型本身的问题。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询