LTI系统时域分析:Matlab建模、卷积与稳定性验证

发布时间:2026/9/17 16:44:47
LTI系统时域分析:Matlab建模、卷积与稳定性验证 简介本资源是华南理工大学《信号与系统》课程配套的第二份实验报告面向电子信息、自动化等专业本科生聚焦离散信号频谱分析的核心能力训练。报告系统讲解DFT/FFT原理、参数选取如点数N对分辨率与泄漏的影响、误差来源识别混叠、截断、采样率不足及改善策略补零、非矩形窗设计、周期整数倍截取并通过8组MATLAB实操任务展开验证涵盖cosine序列、指数衰减信号、多频叠加信号等典型场景并附完整代码、频谱/相位图及对比分析。资源为1个480KB的docx文档结构清晰含实验目的、内容、详细代码实现、图像结果与深度讨论便于课后复盘、课程设计参考或考研复试准备。目前已有302人学习下载是理解数字信号处理中频域分析方法的高质量教学实践材料。1. 信号与系统实验报告二不是交作业的模板而是理解LTI系统行为的实操入口在华南理工大学《信号与系统》课程中“实验报告二”通常对应连续时间LTI系统的时域分析核心环节——它不只考察Matlab绘图能力更检验你能否通过冲激响应、阶跃响应、卷积运算和微分方程求解四条路径交叉验证同一系统特性。很多学生把报告写成命令堆砌却没意识到当impulse(sys)和conv(h,x)给出不一致结果时问题不在代码语法而在采样间隔Ts是否满足奈奎斯特准则、初始条件是否被隐式忽略、或连续系统离散化时采用的是零阶保持还是双线性变换。本报告真正价值在于建立“数学模型→数值实现→物理可测性”的闭环思维。适合已完成拉普拉斯变换推导、能手算一阶/二阶系统响应、但尚未打通仿真与理论映射关系的本科生对研究生而言它是重构经典控制建模直觉的最小可靠基线。2. 用Matlab构建连续LTI系统模型从微分方程到状态空间的三重等价验证2.1 微分方程直接建模tf()函数的隐含假设必须显式声明实验报告二常见题型是给定二阶微分方程$$\frac{d^2y(t)}{dt^2} 3\frac{dy(t)}{dt} 2y(t) x(t)$$直接调用sys tf([1], [1,3,2])看似正确但该命令默认将系统视为零初始条件下的传递函数模型。若题目要求“初始状态$y(0)1, \dot{y}(0)0$”则tf()无法描述——此时必须转向状态空间模型。关键参数说明分子向量[1]对应输入导数项系数分母[1,3,2]按$s^2,s^1,s^0$降幂排列顺序错误会导致极点符号反转。% 正确构建明确标注采样属性虽为连续系统但后续离散化需此声明 sys_tf tf([1], [1,3,2], InputName,u,OutputName,y); % 验证查看极点确认稳定性 pole(sys_tf) % 输出应为[-1,-2]实部全负→渐近稳定提示tf()创建的对象默认Ts0连续时间但若后续误用c2d(sys_tf,0.1)进行离散化未指定方法时默认采用零阶保持ZOH而ZOH在高频段会引入相位滞后——这正是实验中阶跃响应超调量偏大的常见原因。2.2 状态空间建模显式承载初始条件的物理建模方式将微分方程转化为状态变量形式令$x_1y, x_2\dot{y}$得$$\begin{cases} \dot{x}_1 x_2 \ \dot{x}_2 -2x_1 -3x_2 u \end{cases}$$对应矩阵$$A \begin{bmatrix}01\-2-3\end{bmatrix},\ B \begin{bmatrix}0\1\end{bmatrix},\ C \begin{bmatrix}10\end{bmatrix},\ D 0$$A [0,1; -2,-3]; B [0;1]; C [1,0]; D 0; sys_ss ss(A,B,C,D); % 设置非零初始状态实验报告二常考场景 x0 [1;0]; % y(0)1, dy/dt(0)0 t 0:0.01:5; u ones(size(t)); % 单位阶跃输入 [y_ss,tout,x] lsim(sys_ss,u,t,x0); % 注意lsim自动处理初始状态2.2.1 三重模型一致性验证为什么impulse()和lsim()结果必须重合理论上冲激响应$h(t)$是系统对$\delta(t)$的零状态响应而lsim()在x0[0;0]且udirac(t)时应与impulse(sys)完全一致。但实际中dirac(t)在离散时间无法表示Matlab用u(1)1/Ts, u(2:end)0近似Ts为采样间隔若Ts0.01则近似冲激强度为100导致lsim()输出幅值放大100倍解决方案统一用impulse()生成基准再用lsim()验证零状态响应% 获取精确冲激响应连续时间解析解 [h_imp,t_imp] impulse(sys_ss,t); % 用lsim模拟零状态冲激响应输入设为窄脉冲 u_imp zeros(size(t)); u_imp(1) 1/0.01; % 保证∫u dt ≈1 [y_lsim,t_lsim] lsim(sys_ss,u_imp,t,[0;0]); % 绘图对比需归一化 plot(t_imp,h_imp,b,t_lsim,y_lsim,r--); legend(impulse(),lsim() with narrow pulse);2.3 传递函数与状态空间的自动转换ss2tf()的精度陷阱当从状态空间转回传递函数时[num,den] ss2tf(A,B,C,D)可能因浮点误差产生虚假高次项[num,den] ss2tf(A,B,C,D); % 可能返回 num[1.0000e00, -1.1102e-16], den[1.0000, 3.0000, 2.0000] % 多余的-1.11e-16项会导致roots(den)计算出虚部1e-8的伪复根 % 正确做法手动清理数值噪声 num round(num*1e10)/1e10; % 保留10位小数精度 den round(den*1e10)/1e10;3. 卷积运算的实操落地离散近似连续卷积的三个关键控制参数3.1conv()函数的底层逻辑为什么必须补零且调整时间轴连续卷积定义为$y(t)\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$离散化后为$y[n]\sum_k x[k]h[n-k]$。但Matlab的conv(x,h)直接返回长度为length(x)length(h)-1的序列其时间轴需手动校准% 假设x(t)为[0,2]上的矩形脉冲h(t)为系统冲激响应 Ts 0.01; % 采样间隔决定频谱混叠风险 t_x 0:Ts:2; x (t_x0 t_x2); % 矩形脉冲 t_h 0:Ts:5; h impulse(sys_ss,t_h); % 获取冲激响应 % 直接conv结果的时间向量 y_conv conv(x,h)*Ts; % *Ts补偿离散积分的Δt t_y (0:length(y_conv)-1)*Ts; % 错误未对齐原点 % 正确时间轴卷积起点为t_x(1)t_h(1)0终点为t_x(end)t_h(end) t_y_correct 0:Ts:(t_x(end)t_h(end)); % 但y_conv长度比t_y_correct多1需截断 y_conv y_conv(1:length(t_y_correct));3.1.1 采样间隔Ts的双重影响精度与计算量的平衡Ts取值时域精度频域混叠风险conv()计算量实验报告常见问题0.001高低带宽500Hz极大内存溢出运行超时、OOM0.02中中若h(t)含25Hz分量可接受阶跃响应振荡失真0.1低高丢失高频细节极小冲激响应峰值偏低30%推荐策略先用bandwidth(sys_ss)获取系统带宽BW再按Ts 1/(10*BW)设定10倍过采样。对本例二阶系统极点-1,-2BW≈2.8 rad/s≈0.45 Hz故Ts≤0.22即可取0.05兼顾精度与效率。3.2 解析解与数值解的交叉验证用dsolve()校准卷积结果对简单输入如单位阶跃可手算或调用符号工具箱求解析解syms t y(t) Dy diff(y); D2y diff(y,2); eqn D2y 3*Dy 2*y heaviside(t); % heaviside即u(t) cond [y(0)0, Dy(0)0]; % 零初始条件 y_analytic dsolve(eqn,cond); y_analytic simplify(y_analytic); % 得到1 - exp(-t) 0.5*exp(-2*t) % 数值化对比 t_num 0:0.05:5; y_num 1 - exp(-t_num) 0.5*exp(-2*t_num); y_conv_interp interp1(t_y_correct, y_conv, t_num, pchip); % 用三次样条插值对齐 max(abs(y_num - y_conv_interp)) % 应1e-3注意heaviside(t)在t0处定义为0.5而数值阶跃uones(1,N)在t0处为1此差异会导致t0点误差。实验报告中若要求严格匹配需在数值输入中设置u(1)0.5。4. 实验报告二的三大高频扣分点及规避方案4.1 冲激响应绘图中的坐标轴陷阱impulse()默认显示范围失效Matlab的impulse(sys)自动选择时间范围但对慢衰减系统如极点接近虚轴默认范围可能截断尾部。例如极点为-0.1±j3时impulse()仅显示前20秒而实际响应需50秒才衰减至1%以下。强制指定时间范围t_custom 0:0.1:100; % 覆盖5τ50秒再加裕量 h_custom impulse(sys_ss,t_custom); plot(t_custom,h_custom); xlabel(t (s)); ylabel(h(t)); % 关键添加理论衰减包络线验证数值精度 tau 1/0.1; % 时间常数 envelope exp(-t_custom/tau); hold on; plot(t_custom,envelope,k--,LineWidth,1.2); legend(h(t),e^{-t/\tau});4.2 卷积结果的归一化争议能量守恒检验不可省略离散卷积yconv(x,h)*Ts满足$\sum y[n]\cdot Ts \approx \int y(t)dt$但学生常忽略*Ts导致能量偏差。验证方法% 输入能量矩形脉冲 E_x sum(x.^2)*Ts; % 离散能量定义 % 系统能量增益L2范数 G norm(h)^2 * Ts^2; % 连续系统能量增益为∫|h(t)|²dt % 输出能量 E_y sum(y_conv.^2)*Ts; % 检查是否满足E_y ≈ E_x * G帕塞瓦尔定理 abs(E_y - E_x * G)/E_y % 应5%4.3 报告结论部分的致命漏洞混淆“系统稳定”与“响应有界”许多报告写道“因冲激响应衰减系统稳定”。这是不严谨的——BIBO稳定要求∫|h(t)|dt∞而衰减只是必要非充分条件。对本例系统% 数值验证BIBO稳定性 L1_norm sum(abs(h_custom))*Ts; % ∫|h(t)|dt的离散近似 if L1_norm Inf L1_norm 0 fprintf(BIBO稳定L1范数%.4f\n, L1_norm); else fprintf(BIBO不稳定\n); end % 对于本例输出应为L1_norm≈1.0有限正值5. 进阶技巧用lsim()反推系统参数——实验报告二的隐藏考点当实验提供实测输入输出数据如示波器截图需从u(t),y(t)反推系统模型。这不是拟合而是利用LTI系统特性5.1 从阶跃响应提取时间常数与阻尼比对二阶系统阶跃响应的超调量Mp和峰值时间tp直接关联阻尼比ζ和自然频率ω_n$$M_p e^{-\pi\zeta/\sqrt{1-\zeta^2}},\quad t_p \frac{\pi}{\omega_n\sqrt{1-\zeta^2}}$$Matlab实现% 假设已获得阶跃响应y_step零初始条件 [y_max, tp_idx] max(y_step); Mp (y_max - 1)/1; % 相对于稳态值1的超调 zeta_est -log(Mp)/sqrt(pi^2 log(Mp)^2); % 阻尼比估计 tp t(tp_idx); wn_est pi/(tp*sqrt(1-zeta_est^2)); % 自然频率估计 % 构建估计模型 sys_est tf([wn_est^2], [1, 2*zeta_est*wn_est, wn_est^2]); % 对比原始响应 figure; step(sys_ss,b,sys_est,r--,5); legend(原系统,估计模型);5.2 利用卷积逆运算识别未知系统若已知输入x(t)和输出y(t)可通过频域除法求H(jω)% 确保x,y等长且去均值 x x - mean(x); y y - mean(y); N length(x); X fft(x); Y fft(y); H_est Y ./ X; % 频域传递函数估计 % 处理X中零点避免除零 H_est(abs(X)1e-10) 0; h_est ifft(H_est); % 截取物理相关部分因果系统h(t)0 for t0 h_est real(h_est); h_est(1:floor(N/2)) 0; % 强制因果性此方法在实验报告二中常作为拓展题给定一组u(t),y(t)测量数据要求写出h(t)表达式并验证。关键在于频域除法前必须加窗如汉宁窗抑制频谱泄漏否则H_est在高频段噪声极大。提示实际测量中y(t)含噪声直接Y./X会放大噪声。工业标准做法是使用tfestimate(u,y)获取互功率谱再计算H(f)S_{yu}(f)/S_{uu}(f)该函数自动应用平均和窗函数——这才是真实系统辨识的起点。本文还有配套的精品资源点击获取