简介：面向控制系统与Simulink入门学习者的MATLAB小工具包，聚焦于把S平面上的传递函数转换为直观的三维曲面图，便于观察频率响应随频率与相位变化的特性，辅助理解极点零点分布、稳定性与动态性能等抽象概念。压缩包共2个文件，以m脚本和txt文本为主，脚本承担传递函数定义、频率范围设置、响应计算与三维绘图等核心逻辑，文本文件为使用许可说明；整体约2KB，轻量易取用。已有1822人学习，说明该主题在课程实验与自学中具有一定关注度。读者可获得一份可直接参考的传递函数三维可视化实现，用于对照tf函数建模、频响计算及surf、meshgrid绘图思路，也可作为Simulink基础学习的补充素材，适合需要把控制理论曲线与三维图形对应起来的中初级学习者。

1. 三张图看清一个传递函数：从 tf 对象到频域曲线

某开发者把 Simulink 里的控制器参数改完，打开 Bode 图一看，幅频曲线往右挪了，相位裕度却从 45 度掉到 12 度，跟手算对不上。查了半天，问题不在模型——他画图时手敲的分母系数漏了一项。

在 MATLAB 里做传递函数绘图，这类翻车太常见：tf、zpk、ss 三种写法混着用，Bode、Nyquist、Nichols 三张图各画各的，中间没有一次交叉验证。

这份资源围绕“一个传递函数、三张频域图”搭了一套 MATLAB 脚本与 Simulink 示例模型，从建模、出图、算稳定裕度，一路走到与 Simulink 线性化结果核对。刚接触 Simulink 基础、需要把频域指标落成具体数字的人能照着跑；已经会画图、被单位、时延、高阶数值问题坑过的熟手，也能在里面找到边界条件。

2. 传递函数的三种写法与互转：tf、zpk、ss 该用哪个

传递函数在 MATLAB 里不是一段公式字符串，而是一个带属性的对象。你敲tf([1 3],[1 4 5 0])，存下来的是分子分母系数、输入输出名、时延、采样时间这一堆字段。后面画 Bode、算裕度、跟 Simulink 做线性化核对，读的全是这些字段。所以第一步不是急着调bode，而是先想清楚用哪种表示把这个对象建出来。

选错了表示，图能画出来，但数字不一定可信。我见过把八阶系统的分母手写展开再喂给 tf 的，图看着挺顺，裕度算出来是负的，仿真却稳得很——这种“图对、数错”的情况，比直接报错难查得多。

2.1 三种写法不是三个东西：选型的三个判断点

tf（多项式）、zpk（零极点增益）、ss（状态空间）描述的是同一个线性时不变系统，差别在数值稳定性和信息密度，不在功能。

低阶、系数整齐、要给别人看分母长什么样，用 tf 最直接。二阶以内的系统手写 tf 几乎不出错，报告里贴出来也一目了然，评审一眼能对上物理量。

阶次过六，或者零极点分布跨了好几个数量级，我会强制换成 zpk。原因很实在：tf 内部要靠多项式乘法展开分母，展开本质是卷积，系数一大一小挨着，浮点有效位掉得很快。zpk 从头到尾只存零极点位置和增益，不展开，救的就是这一段精度。

MIMO 系统、后面要接状态初值、要做状态反馈观测器，用 ss。Simulink 的线性化输出天生就是 A、B、C、D 四个矩阵，转成 tf 属于二次加工，多一道工序就多一次失真机会。

判断顺序我一般这么走：先看输入输出是不是多于一对，是就走 ss；再看阶次是不是超过六，是就走 zpk；剩下的小系统才用 tf。这个顺序不是死规矩，但能挡掉大部分精度问题。

2.2 建立与互转：一段能跑通的脚本

包里第一节的脚本就是干这件事的：建一个三阶系统，派生三种表示，再用频响数据验证互转有没有丢东西。别用肉眼判断“两条曲线看起来一样”，要落到数字上。

% 同一个三阶系统，三种写法 num = [1 3]; % 分子：s + 3 den = [1 4 5 0]; % 分母：s^3 + 4s^2 + 5s，含一个积分环节 G_tf = tf(num, den); % 多项式形式 G_zpk = zpk(G_tf); % 零极点增益形式，由 tf 直接派生 G_ss = ss(G_tf); % 状态空间形式，默认可控标准型 % 互转回来后逐频点比对，不要靠眼睛 w = logspace(-2, 2, 300); H_tf = squeeze(freqresp(G_tf, w)); H_zpk = squeeze(freqresp(tf(G_zpk), w)); H_ss = squeeze(freqresp(tf(G_ss), w)); % 逐点相对误差，用最大相对偏差衡量回转是否无损 err_zpk = max(abs(H_tf - H_zpk) ./ max(abs(H_tf), eps)); err_ss = max(abs(H_tf - H_ss) ./ max(abs(H_tf), eps)); fprintf('zpk 回转误差 %.3e，ss 回转误差 %.3e\n', err_zpk, err_ss);

逻辑上先建 tf 作为基准，再由它派生 zpk 和 ss，保证三者描述的是同一个系统；回转比对用freqresp在频域取复数响应，逐点做相对误差。这里特意没用norm：分母里有 s 因子，系统在原点是极点，H2 范数不存在，直接norm(G_tf - tf(G_zpk))会抛错，新手往往卡在这一步。

参数上，logspace(-2,2,300)覆盖 0.01 到 100 rad/s，对一个三阶系统足够看清全部转折；eps是防止某个频点响应恰好为零时除零；freqresp返回的是 1×1×N 的三维数组，squeeze压成一维才能跟频率向量对齐。误差量级正常在 1e-14 附近，如果跑出来是 1e-3 甚至更大，说明你的分母系数本身就有问题，别急着往下画图。

2.3 离散化与时延：画图之前先变成有理函数

连续系统的图好画，麻烦的是离散控制器和纯滞后。这两类东西直接往 tf 里塞，图能出来，但读出来的裕度不可信。

离散化先要定方法。c2d不写方法时默认零阶保持器，如果 Simulink 模型里的离散模块用的是双线性，两边模型根本不是同一个，后面核对必然对不上。

Ts = 0.01; % 采样周期，必须与 Simulink 里的 FixedStep 一致 Gd_zoh = c2d(G_tf, Ts, 'zoh'); % 零阶保持器，最贴近实际数字控制器 Gd_tus = c2d(G_tf, Ts, 'tustin'); % 双线性变换，低频段吻合好，高频段会翘 Gd_mp = c2d(G_tf, Ts, 'matched'); % 零极点匹配，只在通带内可信 % 时延：不要直接把 exp(-sT) 塞进传递函数 G_delay = G_tf; G_delay.InputDelay = 0.05; % 50 ms 纯滞后 G_pade = pade(G_delay, 4); % 4 阶 Pade 有理近似，才能画 Bode 和 Nyquist

Ts这个参数不是随便填的，它决定了离散后频率轴能画到哪里：上限是奈奎斯特频率pi/Ts，图超过这条线就没有物理意义。pade的第二参数是近似阶数，常用 3 到 5 阶，阶数越高在低频段越准，但会引入额外的右半平面零点，高频段相位反而失真更厉害。判稳时只信 Pade 近似带宽以内的那一段，超出的部分画出来好看，不能拿来下结论。

2.4 阶次上去以后：poly 展开的精度账

三个一阶环节串联，手写分母看起来风平浪静，实际系数动态范围能跨三个数量级。阶次再往上走，问题就不是“图丑不丑”，而是“分母还代不代表原系统”。

s = tf('s'); G1 = 1/(s + 1); G2 = 1/(s + 10); G3 = 1/(s + 100); G_conv = G1 * G2 * G3; % 内部按多项式卷积展开分母 G_zpkm = zpk(G1) * zpk(G2) * zpk(G3); % 零极点直接相乘，不展开 [num_c, den_c] = tfdata(G_conv, 'v'); den_nz = den_c(den_c ~= 0); spread = max(abs(den_nz)) / min(abs(den_nz)); % 分母系数的动态范围 fprintf('分母系数动态范围 %.2e\n', spread);

这段代码的重点是最后那个spread。它超过 1e6 就说明分母系数在双精度下已经开始互相“吃掉”有效位，此时pole(G_conv)求出来的极点位置会有可见偏差，再拿去画 Bode、算裕度，误差会一路放大。tfdata(..., 'v')里的'v'是让返回值直接给数值向量而不是元胞数组，省一次cell2mat。真遇到动态范围大的系统，做法是全程用 zpk 建模、用 zpk 相乘，只在导出报告给人看的时候临时转一次 tf。

3. Bode 图落地：频率范围、幅值单位与稳定裕度怎么设

Bode 图是这三张图里最容易被当成“随便一画”的一张。默认频率范围常常把关键转折点挤在最左边一小段，幅值单位默认是 dB 但有人拿线性值去读，相位默认折回 ±180 度导致累计相位看不见。这三件事任意一件出错，读出来的裕度就是错的。

包里第三节的脚本把出图拆成四步：先算频率范围，再设绘图选项，然后算裕度，最后做多曲线对比和数据导出。顺序不能反，频率范围定错了，后面所有读数都建立在错误的坐标系上。

3.1 频率范围不是拍脑袋：由极点模值定上下限

默认的bode(G)会自动选范围，但对含积分环节的系统，它经常从 1e-2 起步，把低频段的 -20 dB/dec 斜率截掉一大块。手算过一次带宽之后，我基本不再用默认范围。

p = pole(G_tf); p = p(abs(p) > 1e-9); % 剔除原点上的积分极点，否则 min 会取到 0 if isempty(p) w_lo = 1e-3; % 纯积分系统，下限手工给 else w_lo = 0.1 * min(abs(p)); % 低于最慢极点一个数量级 end w_hi = 10 * max(abs(p)); % 高于最快极点一个数量级 w = logspace(log10(w_lo), log10(w_hi), 500);

逻辑是让频率窗口把全部极点模值包在中间，两边各留一个数量级。低于最慢极点一个数量级，能看到低频渐近线；高于最快极点一个数量级，能看到高频滚降的斜率对不对。abs(p) > 1e-9这个过滤很关键，积分环节的极点数值上是 0 或者极小的负数，不剔掉会让log10(w_lo)直接变成负无穷。

logspace的第三个参数是点数，500 点对一张单系统 Bode 图足够密，做数据导出时我会加到 2000 点，避免 CSV 里曲线出现折线感。

3.2 bodeoptions 里真正影响结论的四个参数

bodeoptions有一堆属性，绝大多数不影响你读出来的数字，但这四个会。

opts = bodeoptions; opts.FreqUnits = 'rad/s'; % 与推导时用的单位一致，别在 Hz 和 rad/s 之间来回换 opts.MagUnits = 'dB'; % 读裕度必须用 dB，线性刻度容易看错量级 opts.PhaseUnits = 'deg'; opts.Grid = 'on'; opts.XLimMode = 'manual'; opts.XLim = [w_lo w_hi]; opts.PhaseWrapping = 'off'; % 不折回 ±180，看清累计相位 bodeplot(G_tf, opts);

FreqUnits是最容易出事的一个。Simulink 里很多模块的频率参数按 Hz 填，脚本里按 rad/s 推，混一次就差 2π 倍，穿越频率能差出一个数量级。PhaseWrapping设成'off'之后相位会一路往下走，三阶系统能看到 -270 度，这在判断“相位还差多少到 -180”时比折叠显示直观得多。

XLimMode必须显式设成'manual'，只给XLim而不改模式，MATLAB 会按自动范围画，你给的区间被忽略。这个坑我踩过一次，调了半天以为选项名字写错了。

3.3 用 margin 把曲线读成数字

图是给人看的，数字是给判据用的。margin一次给出四个量：幅值裕度、相位裕度、相位穿越频率、幅值穿越频率。

[Gm, Pm, Wcg, Wcp] = margin(G_tf); Gm_dB = 20*log10(Gm); % margin 返回的是线性幅值比，报告里写 dB if isinf(Gm) fprintf('幅频不穿越 0 dB，幅值裕度无穷大，Wcp = %.3f rad/s\n', Wcp); end T = feedback(G_tf, 1); % 单位负反馈闭环 bw = bandwidth(T); % 闭环 -3 dB 带宽 fprintf('Gm = %.2f dB, Pm = %.2f deg, Wcg = %.3f, Wcp = %.3f, BW = %.3f\n', ... Gm_dB, Pm, Wcg, Wcp, bw);

margin返回的Gm是线性比值不是 dB，这一步换算漏掉的话，报告里的幅值裕度会写成 3.16 而不是 10 dB，差了一个对数关系。Wcg在一阶系统和过阻尼系统上常常是Inf或者NaN，因为相位根本到不了 -180 度，此时幅值裕度是无穷大，属于正常现象，不是脚本报错。

bandwidth要喂闭环，不是开环。开环的bandwidth在 MATLAB 里没有物理意义，喂错了会得到一个跟闭环带宽差很远的数。

3.4 多条曲线同图对比与数据导出

调参阶段真正有用的不是单条曲线，是同一坐标系下三四个增益版本的对比。

G_list = [G_tf, G_tf*1.2, G_tf*0.8]; % 标称、增益 +20%、增益 -20% h = bodeplot(G_list, opts); legend('标称', '增益 +20%', '增益 -20%', 'Location', 'southwest'); % 把幅频和相频导成表格，方便贴进报告或做后处理 [mag, phase, wout] = bode(G_tf, w); T_out = table(wout(:), 20*log10(squeeze(mag)), squeeze(phase), ... 'VariableNames', {'w_rad_s', 'mag_dB', 'phase_deg'}); writetable(T_out, 'bode_export.csv');

bode三输出的顺序是幅值、相位、频率，跟bodeplot的句柄用法不一样，别混。mag是线性幅值，导数时手动20*log10一下，导出的 CSV 才和图上看到的一致。squeeze在这里不能省，1×1×500 的数组直接进 table 会变成一列元胞。

4. Nyquist 与 Nichols：包围圈判稳和等 M 圆读图

Bode 图告诉你“离 -1 点还有多远”，Nyquist 图告诉你“绕没绕过去”，Nichols 图告诉你“闭环峰有多高”。三张图回答的是三个不同问题，只画一张就下结论，很容易在开环不稳定系统上翻车。

包里第四节的脚本把三张图串起来跑，同一组频率点、同一套单位，画完互相核对。这一节的重点不是画图命令本身，是判据怎么用、参数怎么设才能看清 -1 点附近那一小片区域。

4.1 Nyquist 判稳：Z = P + N，符号和方向最容易记反

Nyquist 判据的完整形式是Z = P + N：Z 是闭环右半平面极点数，P 是开环右半平面极点数，N 是 Nyquist 曲线顺时针包围 -1 点的圈数。闭环稳定等价于 Z 等于零。

麻烦在于 N 的符号约定。逆时针包围算正、顺时针算负，不同教材写法不一样，记串一次结论就反了。我的做法是不背符号，直接用allmargin的Stable字段交叉确认。

P = sum(real(pole(G_tf)) > 0); % 开环右半平面极点数 S = allmargin(G_tf); % 结构体：GainMargin、PhaseMargin、GMFrequency 等 fprintf('开环不稳定极点 %d 个，各裕度对应闭环稳定标志：\n', P); disp(S.Stable); % 手工取一段频率的实部虚部，用于在图上定位 -1 点附近 [re, im] = nyquist(G_tf, w); re = squeeze(re); im = squeeze(im); idx = find(abs(re + 1) < 0.2 & abs(im) < 0.5); % 靠近 (-1,0) 的频点 fprintf('靠近 -1 点的频段：%.3f ~ %.3f rad/s\n', min(w(idx)), max(w(idx)));

allmargin的Stable是逻辑数组，每个元素对应一组裕度下的闭环稳定性判断，全为真才算稳。比手数圈数可靠得多，尤其是曲线在 -1 点附近来回穿插的时候。最后那段找idx的代码是为了定位“该放大哪一段”，下一小节会用到。

4.2 nyquistoptions 与 -1 点附近的放大

默认的 Nyquist 图会把正负频率两支都画出来，镜像分支经常从 -1 点前面横穿过去，把关键区域挡得严严实实。

opts_n = nyquistoptions; opts_n.FreqUnits = 'rad/s'; opts_n.ShowFullContour = 'off'; % 只画正频率段，避免镜像分支遮挡 opts_n.XLimMode = 'manual'; opts_n.XLim = [-2 1]; opts_n.YLimMode = 'manual'; opts_n.YLim = [-1.5 1.5]; nyquistplot(G_tf, opts_n);

ShowFullContour设成'off'只影响显示，不影响判据计算，这一点可以放心。XLim给[-2 1]、YLim给[-1.5 1.5]是把窗口压在 -1 点周围两倍范围，曲线绕没绕过去一眼能看出来。范围开太大，-1 点缩成一个像素，等于白画。

对含积分环节的系统，Nyquist 曲线在低频端会趋向无穷远，此时窗口要给得稍微宽一点，否则曲线从图外飞进来，看不出它从哪个方向绕。

4.3 Nichols 图叠加 ngrid：直接从图上读闭环谐振

Nichols 图的横轴是开环相位、纵轴是开环幅值，叠上ngrid之后，曲线和等 M 圆的切点就是闭环谐振峰，不用再算一遍闭环。

nichols(G_tf); ngrid; % 等 M 圆 + 等相角线 hold on; nichols(G_tf * 1.2); % 增益抬 20% 后的曲线 hold off;

读图方法：曲线与 0 dB 横线的交点对应幅值穿越频率，与 -180 度竖线的距离就是相位裕度，曲线离哪条等 M 圆最近，闭环谐振峰 Mr 就在那一档附近。增益抬高 20% 之后，曲线整体上移，如果它擦到了 3 dB 的等 M 圆，说明闭环超调会明显变大。

ngrid必须在nichols之后调用，顺序反了网格会被曲线盖住。旧版本里ngrid会新开一个图窗，新版直接叠加在当前坐标轴上，跑出来不对先查版本。

4.4 三张图交叉验证：一张对照表

同一组数据画三张图，对不上的地方就是问题所在。包里给了一张对照表，我自己也常用这套判据。

三张图读数如果互相矛盾，九成是单位或者频率范围出了问题。先查FreqUnits，再查MagUnits，最后查是不是把开环和闭环搞混了。这个排查顺序比逐行看代码快得多。

5. 避坑清单：五条把图跑歪的血泪记录

这一节列的五个问题，全部是“脚本不报错、图也画出来了、结论是错的”。它们不会在命令行里给你红色警告，只会安静地把一个错误的裕度写进报告。每条按现象、原因、解决三段记，照着对一遍能省掉大半天。

5.1 幅值单位写成 abs，裕度差了一个对数关系

现象：bodeplot出来的曲线在 0 dB 附近穿过，margin给的Gm却是 3.16，报告里直接抄成“幅值裕度 3.16 dB”。

原因：MagUnits设成'abs'时图按线性刻度画，而margin内部始终按线性幅值比计算，两边一个用线性一个用对数，读图的人把 3.16 当成 dB 了。

解决：绘图选项统一写opts.MagUnits = 'dB'；代码里需要线性值时用10^(Gm_dB/20)换算，别在两种量纲之间来回猜。所有打印语句里带单位的量，变量名后缀都加上_dB。

5.2 时延直接塞进 tf，相位一路卷到 -720 度

现象：给系统加 50 ms 纯滞后，Bode 相位从 -90 度一路掉到 -720 度，Nyquist 图上出现好几圈螺旋，判稳完全没法看。

原因：exp(-sT)不是有理函数，画图时按无穷阶展开处理，相位随频率线性下降，时延越大绕的圈越多。Nyquist 判据里的“包围圈数”被时延产生的额外圈数污染，数出来的 N 是假的。

解决：先pade(G_delay, 4)做有理近似再画图，阶数取 3 到 5；并且只在 Pade 近似有效的带宽内读裕度。判稳时如果时延占比很大，老老实实回到时域仿真看阶跃响应，别硬用频域判据。

5.3 高阶系统用 tf 展开，分母系数已经失真

现象：八阶系统，Bode 图高频段幅值不降反升，margin算出来的相位裕度是负数，Simulink 仿真却稳得很。

原因：poly展开时系数动态范围跨了七八个数量级，双精度有效位不够，分母多项式本身已经不代表原系统，pole求出来的极点位置都是错的。

解决：从建模开始就用 zpk，环节相乘也用 zpk 相乘，全程不展开；只在导出报告给人看的时候临时转一次 tf。判断标准就是第 2 章那个spread，超过 1e6 就别再用 tf 做计算。

5.4 margin 的输入给成了闭环

现象：裕度算出来永远是无穷大，或者相位裕度九十多度，看着特别“稳”。

原因：margin要的是开环传递函数。把feedback(G,1)的结果喂进去，环路已经闭合，不存在“穿越 0 dB”这回事，返回的Inf是数学上的必然，不是系统真的好。

解决：margin(G_open)，G_open是控制器乘被控对象；闭环只用来算bandwidth、看阶跃响应、算灵敏度函数。脚本里给变量起名时就把_open和_cl分开，别用G1、G2这种名字。

5.5 离散系统没指定 c2d 方法，图和仿真对不上

现象：脚本画的 Bode 和 Simulink 里 Linear Analysis 出的曲线，在 0.5 倍奈奎斯特频率附近差了 3 dB 以上，低频段却几乎重合。

原因：c2d不写方法时默认'zoh'，而 Simulink 里的离散模块用的是双线性或者实际控制器就是双线性离散，两边根本不是同一个模型。低频段两种方法吻合，高频段才分家，所以容易漏看。

解决：c2d显式写方法，和 Simulink 模型里的离散模块设置逐项对齐；核对时两边用同一个Ts。图只看奈奎斯特频率以内的部分，超过pi/Ts的曲线是数学延拓，没有物理意义。

6. 从脚本到 Simulink 联调：用 linmod 校核三条曲线

脚本里的传递函数是手写的，Simulink 模型是搭出来的，两者是不是同一个东西，得用数据说话。我一般用linmod在工作点附近把模型线性化成状态空间，再和脚本里的 tf 在同一个频率向量上比频响。

load_system('plant_demo'); % 打开示例模型，不弹窗口 sys_lin = linmod('plant_demo'); % 工作点附近线性化，返回 a/b/c/d 结构体 G_lin = ss(sys_lin.a, sys_lin.b, sys_lin.c, sys_lin.d); w = logspace(-2, 2, 300); H_script = squeeze(freqresp(G_tf, w)); H_lin = squeeze(freqresp(G_lin, w)); dev_dB = 20*log10(abs(H_lin - H_script) ./ max(abs(H_script), eps)); [dev_max, i_max] = max(dev_dB); fprintf('最大偏差 %.2f dB @ %.2f rad/s\n', dev_max, w(i_max));

linmod只对连续模型有效，模型里只要有离散模块就会报错或者给出错误结果，这种情况要换成linearize配合linio指定分析点。偏差在 -40 dB 以下说明两边模型一致；偏差集中在某个频段，多半是某个环节的饱和或者限幅在工作点附近被激活了，线性化结果不再代表原系统。

调参阶段我还会把裕度扫成一张表：增益从 0.5 倍扫到 2 倍，每步算一次margin，把Pm和bw并排列出来，比一张一张图翻快得多。

从那以后我每次出频域图之前，都强制走一遍三张图核对：Bode 看形状和裕度，Nyquist 看包围方向，Nichols 看闭环峰，三个对不上就先查单位和开环闭环。这套习惯救过我至少两次，一次是分母漏项，一次是c2d方法没对齐。希望帮到你。

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