Matlab连续时间LTI系统实验：时频分析、滤波器设计与音频处理

简介：这份资源是面向通信工程、电子信息工程等专业本科生与研究生的连续时间LTI系统综合实验报告，围绕时域与频域分析、系统函数零极点分布、直接/级联/并联三种实现结构对比以及LTI滤波器设计展开，适合作为课程实验参考或自主学习材料。压缩包内共1个docx文档，约2.62MB，完整记录了实验目的、原理推导、Matlab代码实现与结果分析，涵盖冲激响应、阶跃响应、幅频相频特性、参数微变对稳定性的影响以及音频信号频谱分析等环节。已有1303人学习，读者可据此对照复现仿真流程，理解零极点与频率特性、稳定性之间的关系，并掌握不同结构对参数敏感性的差异，从而加深对信号频谱与系统频率特性等核心概念的理解。

1. 从一份实验报告说起：连续时间LTI系统到底怎么用Matlab跑通

很多同学学完信号与系统，公式推得飞起，一打开Matlab就不知道从哪下手。这份资源就是一份完整的连续时间LTI系统综合实验工程，把时域分析、频域分析、零极点分布、滤波器设计、音频信号处理串成一条线，全部用Matlab脚本落地。它适合通信、电子信息类专业的本科生拿来复现实验，也适合已经工作、想把信号处理基础重新捡起来的开发者当练手项目。整份材料不是零散代码片段，而是从系统函数设计到频响验证、再到真实音频滤波的完整链路，拿到手改几个参数就能跑出自己的结果。

2. 工程结构与技术栈：这份实验包由哪些脚本和模块组成

2.1 目录组织与文件职责

这份资源本质是一份实验报告加配套Matlab脚本，按实验内容可以拆成四个功能块：基础系统分析、参数敏感性对比、周期信号截断频谱分析、音频采集与滤波。实际落地时，我一般会把它整理成下面这样的目录，方便反复调试和交作业：

lti_lab/ ├── basic_system/ │ ├── sys_define.m % 定义H(s)，画h(t)、幅频、相频 │ ├── step_response.m % 阶跃响应与频谱变化 │ └── struct_compare.m % 直接/级联/并联结构对比 ├── spectrum/ │ └── trunc_spectrum.m % 周期信号不同截断长度频谱 ├── filter_design/ │ ├── filter_freq.m % 滤波器频响分析 │ └── filter_conv.m % 卷积求输出y(t)与Y(jw) └── audio/ └── audio_process.m % 录音、读写、FFT频谱

每个脚本独立可运行，变量命名统一用Hs、b、a、w这类信号与系统里常见的符号，读起来和课本对得上。技术栈就是纯Matlab，依赖Symbolic Math Toolbox（符号运算求拉普拉斯逆变换、傅里叶变换）和Signal Processing Toolbox（freqs、tf2sos、zp2tf、pzmap、step）。没有额外第三方库，Matlab R2016b以上基本都能跑，新版把wavrecord换成audiorecorder即可。

2.2 系统函数与零极点设计思路

基础部分的核心是设计一个4到6阶的稳定连续LTI系统，要求至少有一对靠近虚轴的共轭复数极点。资源里给出的系统函数是：

$$H(s)=\frac{2.8s+1.99}{(s^2+0.2s+0.01)(s+1)(s+2)}$$

分母展开后是4阶，极点分别来自二次因子 $s^2+0.2s+0.01$ 和两个实极点 $-1$、$-2$。二次因子的极点是 $-0.1\pm j0$ 附近的一对共轭复根，实部为负、非常靠近虚轴，这正是用来观察“极点靠近虚轴导致某频率处幅频出现峰值”的典型布置。所有极点实部为负，系统稳定，符合实验对稳定性的要求。

在Matlab里定义这个系统，用tf或者直接写系数向量都行。我习惯先把分子分母系数算出来，后面freqs、step、pzmap全都复用这两个向量：

% 系统函数 H(s) = (2.8s+1.99) / ((s^2+0.2s+0.01)(s+1)(s+2)) % 分母展开：s^4 + 3.2s^3 + 3.41s^2 + 0.62s + 0.02 b = [0 0 0 2.8 1.99]; % 分子系数，高位补零对齐阶数 a = [1 3.2 3.41 0.62 0.02]; % 分母系数，4阶 sys = tf(b, a); % 建立传递函数对象 pzmap(sys); % 画零极点图，确认极点都在左半平面

这里b前面补三个零是因为分子实际只有2阶，而分母是4阶，tf要求两个向量按最高阶对齐。a的常数项是0.02，来自0.01乘1乘2，写代码时容易手算错，建议用conv连乘二次因子和两个一阶因子再取系数，避免翻车：

a = conv([1 0.2 0.01], conv([1 1], [1 2])); % 自动展开分母

这样改极点位置时只改二次因子系数，分母自动重算，比手写展开式可靠得多。

2.3 时域与频域响应的仿真入口

时域响应用符号运算求冲激响应最直观，资源里用ilaplace对Hs做拉普拉斯逆变换，再乘heaviside(t)保证因果性。频域则用fourier求h(t)的傅里叶变换，取模得幅频、取相角得相频。这套写法适合教学演示，因为符号结果和课本公式一一对应；但如果要做大量数值扫描，符号运算会慢，实际工程里更常用freqs直接算频响、用impulse/step算时域。

syms s t w; Hs = (2.8*s + 1.99) / ((s^2 + 0.2*s + 0.01)*(s + 1)*(s + 2)); h = ilaplace(Hs) * heaviside(t); % 冲激响应，乘阶跃保证t<0为0 Hw = fourier(h); % 频域响应 subplot(2,2,1); ezplot(h); title('h(t)'); subplot(2,2,2); ezplot(abs(Hw)); title('|H(jw)|'); subplot(2,2,3); ezplot(angle(Hw)); title('相频');

ilaplace返回的是符号表达式，heaviside(t)把非因果部分压掉，这一步不能省，否则画出来的h(t)在负时间轴上有值，和因果系统定义矛盾。ezplot在新版Matlab里逐渐被fplot替代，如果报错就换成fplot(h, [0 20])，指定绘图区间更可控。

3. 三种实现结构的对比：直接型、级联型、并联型怎么在代码里切换

3.1 用 tf2sos 和 zp2tf 做结构转换

实验要求分别用直接、级联、并联三种结构实现同一个系统，并比较参数微调后零极点、稳定性、时频特性的变化。Matlab里转换结构靠三个函数：tf2sos把传递函数转成二阶节级联（second-order sections），sos2tf转回来，tf2zp和zp2tf在零极点增益形式和传递函数形式之间互转。资源里的做法是先由b,a得到级联和并联系数，再分别画零极点图和频响。

b = [0 0 0 2.8 1.99]; a = [1 3.2 3.41 0.62 0.02]; [sos, G] = tf2sos(b, a); % 直接型 -> 级联型二阶节 [b2, a2] = sos2tf(sos, G); % 级联型 -> 等效直接型系数 [z, p, k] = tf2zp(b, a); % 直接型 -> 零极点增益 [b3, a3] = zp2tf(z, p, k); % 零极点 -> 并联型等效系数 subplot(2,2,1); pzmap(b, a); title('原始零极点'); subplot(2,2,2); pzmap(b2, a2); title('级联结构零极点'); subplot(2,2,3); pzmap(b3, a3); title('并联结构零极点');

sos是一个矩阵，每行是一个二阶节的分子分母系数，G是总增益。sos2tf把多个二阶节乘回去，理论上应该还原出原来的b,a，数值上会有微小舍入误差，这正是后面观察“结构对参数敏感度”的切入点。zp2tf由零极点重建多项式，并联结构在实现上就是把部分分式展开后的各一阶/二阶项相加，Matlab里用residue更直接，但资源用zp2tf做等效验证也够用。

3.2 参数微调与敏感度观察

实验的关键动作是“微量改变系统中的一个或部分参数”，变化量建议在万分之一到5%之间，然后看三种结构的零极点漂移、稳定性、时频响应和阶跃响应差异。资源里把分母常数项从0.02改成0.01，相当于动了靠近虚轴那对共轭极点的阻尼，观察直接型系数变化对整体极点的牵动。

a1 = [1 3.2 3.41 0.62 0.01]; % 只改常数项，模拟参数漂移 w = -4:0.01:4; h0 = freqs(b, a, w); % 原始频响 h1 = freqs(b, a1, w); % 直接型参数改变后频响 subplot(2,1,1); plot(w, abs(h0), w, abs(h1)); legend('原始','参数改变'); title('幅频对比'); subplot(2,1,2); plot(w, angle(h0), w, angle(h1)); title('相频对比');

freqs的第三个参数是频率向量，单位rad/s，范围取-4:0.01:4是因为系统主要频段集中在低频，极点靠近虚轴导致峰值出现在很低的频率处。直接型系数改动会同时影响所有极点位置，级联型每个二阶节只对应一对零极点，改一个节的系数不会波及其他节，这就是实验想说明的“级联/并联结构独立调整零极点、误差积累小”的结论。实际写代码时，级联型要改某个二阶节的系数，直接改sos矩阵对应行再sos2tf回去，比改展开后的a向量更符合结构本意。

3.3 阶跃响应对比与稳定性判断

阶跃响应用step(b,a,t)一行搞定，时间向量取0:0.01:80是为了让靠近虚轴的慢极点充分衰减、曲线进入稳态。四种情况并排画：原始、直接型改参数、级联型改参数、并联型改参数。

t = 0:0.01:80; subplot(2,2,1); step(b, a, t); title('原始阶跃响应'); subplot(2,2,2); step(b, a1, t); title('直接型参数改变'); subplot(2,2,3); step(b2, a2, t); title('级联型'); subplot(2,2,4); step(b3, a3, t); title('并联型');

判断稳定性看极点实部是否全为负，代码里可以用max(real(roots(a))) < 0做自动检查，比肉眼看零极点图可靠。参数改动后如果某个极点越过虚轴，step曲线会发散，这时候别急着怀疑代码，先roots(a1)看极点实部符号，大概率是参数改大了。实验报告里强调“右半平面有极点则不稳定”，代码验证就落在这行判断上。

4. 周期信号截断频谱与LTI滤波器设计：从理论到卷积验证

4.1 不同截断长度的频谱差异

拓展部分构建 $f(t)=\cos(t)+\cos(2t)$，分别截取-10:0、-5:0、-10:10三个区间做数值傅里叶变换，比较幅度谱和相位谱。核心代码是离散求和近似连续傅里叶变换：

r = 0.02; % 采样间隔 t = -10:r:10; % 截断区间，可换成 -10:0 或 -5:0 N = 200; W = 2*pi; k = -N:N; w = k*W/N; % 频率轴 f = cos(t) + cos(2*t); F = r * f * exp(-1j*t'*w); % 数值积分近似傅里叶变换 subplot(3,1,1); plot(t, f); title('f(t)'); subplot(3,1,2); plot(w, abs(F)); title('幅度谱'); subplot(3,1,3); plot(w, angle(F)*180/pi); title('相位谱(度)');

r是时间步长，F = r * f * exp(-1j*t'*w)这一行是矩形法数值积分，t'转置后与w做矩阵乘法，得到每个频点上的积分值。截断长度越短，频谱泄漏越严重，幅度谱上主瓣变宽、旁瓣抬高，这正是频域卷积定理的体现：时域乘矩形窗等于频域卷积sinc函数。三个区间对比时，-10:10因为包含更多周期，谱线最接近理论冲激；-5:0只覆盖两个多周期，泄漏最明显。相位谱在幅度接近零的频点会跳变，这是angle函数在数值噪声下的正常现象，分析时只看幅度较大的频段。

4.2 零极点布置设计滤波器

滤波器设计部分给出 $H(s)=\frac{s^2-0.2s+1.97}{s^2+0.2s+1.01}$，零点在 $0.1\pm j$ 附近、极点在 $-0.1\pm1.4j$ 附近。零点靠近虚轴会在对应频率处形成幅频谷值，极点靠近虚轴形成峰值，两者配合把某个频率分量压下去、保留另一个。用freqs画频响验证：

b_f = [1 -0.2 1.97]; % 滤波器分子 a_f = [1 0.2 1.01]; % 滤波器分母 w = -6:0.01:6; H = freqs(b_f, a_f, w); subplot(2,1,1); plot(w, abs(H)); title('|H(jw)|'); subplot(2,1,2); plot(w, angle(H)); title('相频');

看幅频曲线在 $\omega=1$ 和 $\omega=2$ 附近的取值，就能判断它滤掉了哪个分量。零点位置决定谷值频率，极点位置决定峰值频率，调b_f、a_f的系数本质就是挪零极点。这里有个容易忽略的点：freqs的频率轴单位是rad/s，而cos(t)的角频率是1 rad/s、cos(2t)是2 rad/s，对照时别把Hz和rad/s搞混。

4.3 卷积求输出并验证滤波效果

最后一步用拉普拉斯变换求 $Y(s)=H(s)F(s)$，再逆变换得 $y(t)$，画波形和频谱确认滤波生效。资源里用符号运算：

syms s t; Hs = (s^2 - 0.2*s + 1.97) / (s^2 + 0.2*s + 1.01); f = cos(t) + cos(2*t); Fs = laplace(f); Ys = Hs * Fs; yt = ilaplace(Ys) * heaviside(t); ezplot(yt, [0 6]); title('y(t)');

符号解能直接看出稳态分量里哪个频率被保留，但表达式可能很长，实际调试时我更推荐数值卷积：先impulse或ilaplace取h(t)的数值序列，再和f(t)采样序列做conv，乘r修正积分步长。两种方法结果对得上，就说明滤波器设计没问题。稳态输出取时间轴后半段做FFT，看Y(jw)的峰值位置，如果只剩一个频率分量，滤波目的就达到了。

5. 音频信号实战：录音、频谱分析与常见翻车点排查

5.1 录音与读写的最小可用代码

应用部分把前面学的频谱分析落到真实音频上。新版Matlab用audiorecorder替代老的wavrecord，录3秒单声道、8kHz采样，存成wav再读回来做FFT：

rec = audiorecorder(8000, 16, 1); % 8kHz, 16bit, 单声道 recordblocking(rec, 3); % 录3秒，阻塞式 y = getaudiodata(rec); % 取采样数据，范围[-1,1] audiowrite('luyin.wav', y, 8000); % 写文件 [y, fs] = audioread('luyin.wav'); % 读回 sound(y, fs); % 回放 Y = fft(y, 8000); % 8000点FFT subplot(2,1,1); plot(y); title('波形'); subplot(2,1,2); plot(abs(Y)); title('频谱');

audiorecorder三个参数分别是采样率、位深、通道数，8kHz满足语音带宽又不会让FFT点数太大。fft(y, 8000)里第二个参数指定变换点数，y不足8000点会补零、超过会截断，补零会让频谱看起来更平滑但不会提高真实分辨率，这点在分析时要心里有数。频谱图只画前半段更符合实际，因为实信号FFT对称，后半段是镜像。

5.2 避坑与常见问题排查

现象一：ilaplace报错或返回空。原因多半是Hs里有未定义的符号变量，或者分母多项式没展开成s的幂次形式。解决：先用syms s声明，再用collect或simplify整理Hs，确认分母是多项式而不是因子乘积形式再调ilaplace。

现象二：freqs画出来的幅频和fourier符号结果对不上。原因是freqs的频率轴是rad/s、范围没覆盖到极点峰值，或者符号fourier默认双边谱而freqs只算正频率。解决：把w范围放宽到-10:0.01:10，对比时统一看双边或单边，别拿单边谱和双边谱直接比高低。

现象三：tf2sos报“分母必须为实数”。原因是a向量里混入了符号变量或复数。解决：确保b,a是数值向量，符号推导的结果用double(coeffs(...))转成数值再进tf2sos。

现象四：录音回放没声音或全是噪声。原因是采样率设得和声卡不匹配，或者getaudiodata取到的数据类型不对。解决：采样率用8000或44100这类标准值，取数据后y = y / max(abs(y))归一化到[-1,1]，再sound输出。

现象五：阶跃响应曲线发散，以为系统不稳定。原因是参数微调时把极点推过了虚轴，或者时间向量太长放大了数值误差。解决：roots(a)看极点实部，全负才稳定；时间向量按最慢极点的时间常数取5到10倍即可，不必取到80那么长。

6. 二次开发与改造思路：把实验脚本变成可复用的分析工具

拿到这份资源后，最直接的改造是把散落的脚本收成一个函数库，每个功能一个函数、参数从外部传，避免每次改系统都要动脚本里的硬编码。比如把系统定义抽成make_sys(b,a)返回tf对象，把频响分析抽成plot_freq(sys, wrange)，把结构对比抽成compare_struct(b,a,delta)，主脚本只负责调用和画图。这样换一个系统函数、换一组零极点，只改一行参数就能重跑全套分析。

第二个改造方向是把符号运算和数值运算分开。符号部分只用来推导公式、验证理论结果，数值部分用freqs、impulse、step、fft做批量扫描。我一般会写一个verify_symbolic_vs_numeric的小脚本，同一组参数下两种方法各算一遍，画在同一张图上对比，误差超过阈值就报警。这个习惯在调滤波器参数时特别有用，符号解告诉你理论峰值在哪，数值解告诉你实际代码算出来在哪，两者一对照，参数写错立刻暴露。

第三个方向是音频部分的可视化增强。原始代码只画波形和FFT幅度，实际做二次开发时可以加短时傅里叶变换看频谱随时间变化，或者加一个简单的滑动平均滤波对比滤波前后波形。spectrogram函数一行就能出时频图，配合前面的滤波器设计，能直观看到某个频段被压掉的效果。注意spectrogram的窗长和重叠参数会影响时间分辨率与频率分辨率的权衡，窗越长频率分辨率越高、时间定位越模糊，按信号特点调。

改造时有个血泪经验：别一上来就重构全部脚本，先把能跑通的原版跑一遍、存一份基准输出，再逐块替换成函数，每换一块就和基准对比一次。我见过太多人重构到一半发现某条曲线对不上，回头找是哪次改动引入的，没有基准就只能从头再来。从那以后我每次动这类实验代码，都强制先跑一遍原版存图，改完再逐图比对，确认无误才继续下一步。希望这份拆解能帮你把这份实验资源真正跑起来、改出自己的版本。

本文还有配套的精品资源，点击获取