大家有没有过这种经历：控制理论课上，老师让画系统开环传递函数的对数幅频渐近特性曲线，坐标纸摊开，转折频率一个一个算，从低频段开始斜率一层一层叠，叠到后面怀疑人生。反正我是有的，而且不止一次。后来我在某高校做课程答疑，发现学生最大的问题反而不是“不会算转折频率”，而是转折频率算对了之后，后面整条折线的斜率变化太容易乱。再加上手画半对数坐标纸，刻度小、线多，画完还得擦掉重来。于是我就写了一个MATLAB脚本 bd_asymp.m，专门干一件事：输入一个传递函数，自动把对数幅频渐近特性曲线画出来，每一段的斜率值直接标在线上，转折频率处还带数据点。这篇文章就把脚本的用法、原理、完整代码和踩过的坑一次性说清楚。
1. 为什么需要这个脚本：折线背后的真实需求
1.1 手画幅频渐近线的三个痛点
先说说手画这条曲线到底难在哪。对数幅频渐近特性曲线本质上是一系列分段直线，转折频率由系统的零点、极点位置决定。听起来简单，实际画的时候有三个地方特别容易出问题。
第一是半对数坐标纸。频率轴是对数刻度，要找到某个转折频率在坐标纸上的横向位置，本身就容易看错。频率跨度稍大一点，从0.1到100有三个数量级，贴错一个位置整条线就歪了。
第二是斜率叠加。低频段初始斜率由积分环节个数决定，每经过一个转折频率，斜率要加上或减去20dB/dec或40dB/dec。系统稍微复杂点，零点极点加起来五六个，画到高频段的时候还得回头检查前面叠得对不对。
第三是转折频率重合。比如一个零点和一个极点落在同一个频率上，斜率变化要先合并再叠加，漏掉一个符号整段就反了。这些操作在纸上做起来非常琐碎，而且一旦中间发现错误，前面的图十有八九要作废重画。
1.2 渐近线不是“近似”，而是控制系统设计的骨架
有人可能会说：现在软件都能画精确幅频曲线了，为什么还要跟折线较劲？这个想法我一开始也有，后来教了几年课、做了几个实际项目才反应过来，渐近线反倒是更接近“设计语言”的东西。
做稳定性分析的时候，判断系统稳不稳，看的是穿越频率和相位裕度。穿越频率基本就是幅频曲线与0dB轴的交点，这在折线上读非常方便。做校正装置设计的时候，工程师往往先在Bode图上画出期望的开环幅频特性，再反推控制器结构。期望特性通常直接用渐近线表示，因为精确曲线在转折频率附近的圆滑过渡反而会干扰设计意图。
换句话说，精确曲线用来“验证”性能，渐近线用来“设计”系统。这个脚本画的虽然是折线，但在教学、考试和前期方案阶段，它就是最好用的那套工具。
1.3 为什么不用现成的bode命令
MATLAB自带的bode命令当然能画幅频曲线，而且画得很准。但它画的是精确曲线，不是渐近线。MATLAB长期没有直接输出“渐近折线”的命令，做控制的人要么自己写，要么就靠坐标纸手画。
这个脚本存在的意义，就是把“零极点增益 → 转折频率 → 分段折线 → 斜率标注”这条链路自动打通。你给它一个传递函数，它还你一条带斜率标注的渐近线图，过程完全透明，中间变量都在工作区里，方便二次处理。
2. 脚本用法：从传函到折线只要三行
2.1 输入格式与模型转换
bd_asymp.m 的输入是MATLAB控制系统工具箱里的LTI对象，tf、zpk、ss都行。脚本内部有个关键动作：无论你传进来什么格式，先统一转成zpk模型。
为什么非要zpk？因为画渐近线需要的三样东西——零点、极点、增益——在zpk模型里一目了然。tf格式是多项式比值，要额外解方程求根才能拿到零极点；ss格式更要多绕一步。转换成zpk之后，零极点直接列出来，转折频率就是零极点的模，斜率增量看阶次就能定。
调用方式有三种：
bd_asymp(G); % 自动选择频率范围 bd_asymp(G, 0.01, 100); % 手动指定频率范围 bd_asymp(G, [0.01, 100]); % 兼容写法，向量传范围如果不给频率范围，脚本会按最小转折频率除以10、最大转折频率乘以10来定范围，覆盖所有转折点，高低频各留一个十年。
2.2 完整运行示例
下面这个例子我经常拿来演示，因为它同时包含了一阶环节、积分环节和二阶振荡环节，非常典型。
s = zpk('s'); G = 40 * (s + 2) / (s * (s + 1) * (s + 4) * (s^2 + 2*s + 4)); bd_asymp(G);这个系统的零极点分布是：零点s=-2，极点s=0、s=-1、s=-4，还有一对复极点s=-1±j√3。复极点的自然频率wn=sqrt(1^2+3^2)=2，阻尼比ζ=0.5。
运行之后，脚本会画出从0.1到40左右的对数幅频折线，并且在每一段中点标出当前斜率。对于这个系统，转折频率分别为1、2、4，其中最精彩的是w=2这个点：零点提供+20，复极点对提供-40，合并后是-20。整个图形从低频-20dB/dec起步，依次变为-40、-60、-80，层次非常清楚。
2.3 输出内容怎么读
图上有三样信息：蓝色折线、圆圈数据点、斜率标注文字。圆圈标在转折频率的折点处，方便你对照检查脚本计算的转折频率和手算结果是否一致。斜率文字直接标在每段线中间，读图的时候不用再自己心算。
命令行窗口会额外输出欠阻尼提示。比如上面的例子阻尼比ζ=0.5，离峰值修正的临界值0.707很近，脚本会打印类似“检测到欠阻尼复极点，阻尼比0.50，实际幅频曲线在w=2附近有微小峰值”的提示。这个提示不是画在图上的，因为画在图里会干扰折线的教学效果，打印在命令行里刚好够用。
3. 核心原理拆解：零极点如何变成折线
3.1 环节类型、转折频率与斜率增量对照表
画渐近线的底层逻辑，是把传递函数拆成若干个典型环节的串联，然后利用对数的加法性质，把每个环节的幅频贡献逐段叠加。比例环节只把整条曲线上下平移20lg|K|，不改变斜率；转折频率和斜率变化只由零点、极点对应的环节产生。
| 环节类型 | 对应零极点 | 转折频率 | 斜率增量 |
|---|---|---|---|
| 一阶惯性（极点） | 负实极点 p | wc = | p |
| 一阶微分（零点） | 负实零点 z | wc = | z |
| 二阶振荡（共轭极点） | 复极点 p | wc = | p |
| 二阶微分（共轭零点） | 复零点 z | wc = | z |
| 积分环节 | 原点极点 | 无 | 初始 -20 dB/dec |
| 微分环节 | 原点零点 | 无 | 初始 +20 dB/dec |
每种环节的转折频率，数值上等于零极点距离原点的模。负实极点-1，转折频率就是1；复极点-1±2j，转折频率就是sqrt(1+4)约等于2.236。积分类原点处的零极点不产生转折频率，它负责设定低频初始斜率。
3.2 复共轭对的配对与阻尼比修正
复共轭对是最容易写错的地方。zpkdata返回的极点里，一对共轭复根以两个复数的形式出现，如果脚本不去重，会把同一个转折频率算两次，斜率增量再翻一倍，整条曲线直接歪掉。
我采用的策略是只取“上半平面”的复根，也就是虚部大于0的那些。每一对复共轭根有且只有一个在上半平面，取出来正好每个代表一对，不会重复。这个技巧简洁可靠，实测下来比“两两配对找共轭”要稳得多。
阻尼比修正是另一个值得展开的细节。二阶环节的渐近线在wn处是折点，但实际幅频曲线在wn附近有一个圆滑过渡。对于共轭极点，实际峰值比渐近线高出20lg(2ζ)分贝；对于共轭零点，方向相反。只有ζ小于0.707时这个偏差才明显，这叫欠阻尼状态。脚本会在检测到欠阻尼复极点时打印修正提示，方便你手工在图纸上补一笔。
3.3 相同转折频率的斜率合并
转折频率相同的零点和极点，在图上发生在同一个横坐标位置。比如系统里有个零点s=-2和一个复极点对，而复极点的wn也正好等于2，那么w=2处会出现+20和-40同时存在的局面。
物理上，这两个环节同时在w=2开始改变斜率，折线在同一个x位置只能拐一次，所以必须把斜率增量先合并：+20加上-40等于-20，表示经过w=2后斜率净下降20dB/dec。
脚本里用unique配合round来处理这个问题。先把所有转折频率四舍五入到6位小数，unique自动去重排序，再对每个唯一频率把所有对应的斜率增量求和。这个思路简单，但非常关键，少了这一步同一个转折频率会画出两条竖线段，图形完全没法看。
3.4 初始斜率与频率范围的自动选取
初始斜率由原点处极点数减去原点处零点数决定。每多一个原点极点，低频段就多-20dB/dec；每多一个原点零点，就多+20dB/dec。这个n0值直接决定折线第一段的方向。
比如一个一型系统，例如直流电机转速环开环传递函数，含一个积分环节，n0=1，初始斜率就是-20dB/dec。如果系统含两个积分环节，初始斜率就是-40dB/dec。如果系统是纯微分环节，n0=-1，初始斜率为+20。
频率范围自动选取时，脚本会先找所有转折频率中的最小值和最大值。低频端点取最小值除以10，高频端点取最大值乘以10，确保每个转折点两侧都有足够的直线段展示。如果没有转折频率，比如纯增益系统，默认画0.01到100的范围，此时曲线就是一条水平直线。
4. 完整代码：bd_asymp.m 源码与注释
4.1 脚本源码
这个函数不长，核心逻辑就三段：拆零极点、算转折频率、画折线。直接贴完整代码，复制保存成bd_asymp.m就能用。
function bd_asymp(sys, wmin, wmax) % BD_ASYMP 绘制对数幅频渐近特性曲线（Bode幅频渐近线） % 用法： % bd_asymp(sys) 自动选择频率范围 % bd_asymp(sys, wmin, wmax) 指定频率范围 % bd_asymp(sys, [wmin, wmax]) 指定频率范围（兼容写法） % % 输入 sys 可为 tf、zpk、ss 等 MATLAB 控制系统 LTI 模型。 % 无返回值，在当前图窗绘制折线并标注各段斜率。 % 兼容参数写法 if nargin < 2, wmin = []; end if nargin < 3, wmax = []; end if nargin == 2 && numel(wmin) == 2 tmp = wmin; wmin = tmp(1); wmax = tmp(2); end % 统一转 zpk 模型 if ~isa(sys, 'zpk') try sys = zpk(sys); catch error('输入对象无法转换为 zpk 模型，请检查 sys 是否为合法的 LTI 对象'); end end [z, p, k] = zpkdata(sys, 'v'); z = z(:).'; p = p(:).'; eps0 = 1e-8; % 实零点 / 复零点 z_real = z(abs(imag(z)) <= eps0); z_cplx = z(imag(z) > eps0); % 只取上半平面，每个代表一对共轭零点 % 实极点 / 复极点 p_real = p(abs(imag(p)) <= eps0); p_cplx = p(imag(p) > eps0); % 只取上半平面，每个代表一对共轭极点 % 原点处的零极点 z0 = z_real(abs(real(z_real)) <= eps0); zf = z_real(abs(real(z_real)) > eps0); p0 = p_real(abs(real(p_real)) <= eps0); pf = p_real(abs(real(p_real)) > eps0); % 收集转折频率和斜率增量 wc = []; delta = []; if ~isempty(zf) wc = [wc, abs(real(zf))]; delta = [delta, 20 * ones(1, length(zf))]; end if ~isempty(z_cplx) wc = [wc, abs(z_cplx)]; delta = [delta, 40 * ones(1, length(z_cplx))]; end if ~isempty(pf) wc = [wc, abs(real(pf))]; delta = [delta, -20 * ones(1, length(pf))]; end if ~isempty(p_cplx) wc = [wc, abs(p_cplx)]; delta = [delta, -40 * ones(1, length(p_cplx))]; end % 相同转折频率合并 if ~isempty(wc) [wc, ~, ic] = unique(round(wc, 6)); deltaNew = zeros(size(wc)); for i = 1:length(wc) deltaNew(i) = sum(delta(ic == i)); end delta = deltaNew; else delta = []; end % 初始斜率由原点处极点个数减零点个数决定 n0 = length(p0) - length(z0); slope0 = -20 * n0; % 频率范围 if isempty(wmin) && isempty(wmax) if ~isempty(wc) wmin = min(wc) / 10; wmax = max(wc) * 10; else wmin = 1e-2; wmax = 1e2; end else if isempty(wmin), wmin = 1e-3; end if isempty(wmax) if ~isempty(wc) wmax = max(wc) * 10; else wmax = 1e2; end end end if wmin <= 0 warning('频率下限已自动修正为正数'); wmin = 1e-4; end % 折线计算 freqPts = [wmin, wc, wmax]; Lpts = zeros(size(freqPts)); Lpts(1) = 20 * log10(abs(k)) + slope0 * log10(wmin); curSlope = slope0; for i = 2:length(freqPts) Lpts(i) = Lpts(i-1) + curSlope * (log10(freqPts(i)) - log10(freqPts(i-1))); if i - 1 <= length(delta) curSlope = curSlope + delta(i - 1); end end % 绘图 figure; semilogx(freqPts, Lpts, 'b-o', 'LineWidth', 1.5, 'MarkerSize', 4); grid on; xlabel('频率 \omega (rad/s)'); ylabel('幅值 (dB)'); title('对数幅频渐近特性曲线'); xlim([freqPts(1), freqPts(end)]); % 标注每段斜率 for i = 1:length(freqPts) - 1 if freqPts(i + 1) > freqPts(i) wmid = sqrt(freqPts(i) * freqPts(i + 1)); slopeSeg = slope0; for j = 1:i - 1 slopeSeg = slopeSeg + delta(j); end Lmid = interp1(freqPts, Lpts, wmid); text(wmid, Lmid + 2, sprintf('%d dB/dec', slopeSeg), ... 'FontSize', 9, 'HorizontalAlignment', 'center'); end end % 欠阻尼复极点提示 if ~isempty(p_cplx) zetaList = []; for i = 1:length(p_cplx) wn_i = abs(p_cplx(i)); zeta_i = -real(p_cplx(i)) / wn_i; if abs(zeta_i) < 0.707 && zeta_i > 0 zetaList = [zetaList, zeta_i]; end end if ~isempty(zetaList) zetaStr = sprintf('%.2f, ', zetaList); fprintf('提示：检测到欠阻尼复极点，阻尼比 %s。\n', zetaStr(1:end-2)); fprintf('实际幅频曲线在该转折频率附近存在峰值，修正量约为 +20*lg(2ζ) dB。\n'); end end end4.2 关键设计点说明
代码里有两处设计值得单独说明。第一处是复根取上半平面的技巧，这个在前面原理部分讲过，代码里体现在z_cplx和p_cplx的判断上。这个技巧省掉了繁琐的共轭配对循环，也让后续wm计算天然对。
第二处是折线幅值逐步累加的思路。脚本没有直接调用bode算精确值，而是从低频端开始，用当前斜率乘上对数频率间隔，一段一段往后推。这样画出来的折线严格符合渐近线定义，而且在转折频率处不会出现人为圆角。起点幅值用20lg|K|加上初始斜率与lg(wmin)的乘积，保证了整条曲线的绝对幅值准确。
斜率标注那一段用了两层循环。外层遍历所有线段，内层重新累加斜率。虽然效率不高，但这个脚本画一条系统曲线撑死十几个转折点，完全够用，反而比维护一个斜率数组更直观、不容易出错。
4.3 运行效果与图形解读
拿前面那个示例系统来说，运行后你会看到：低频从w=0.1开始，幅值约20lg40=32dB，然后以-20dB/dec下降；到w=1时，经过极点s=-1，斜率转成-40dB/dec；到w=2时，零点和复极点对同时起作用，净变化-20dB/dec，斜率变成-60；到w=4时通过最后一个极点，斜率变成-80。
整条线在图上呈现为逐段下探的折线，每个转折点都标着圆点，每段中间标着当前斜率数字。对照手算结果，每一段都能对上。这个图形既可以交作业，也可以作为设计阶段的快速估算图。
5. 常见问题与排查技巧实录
5.1 “sys无法转换”或含延迟环节怎么处理
如果输入的不是LTI对象，或者模型里带了时间延迟e^(-τs)，zpk转换会直接报错。时间延迟不是有理传递函数，零极点增益模型表达不了它。
遇到这种情况，我的建议是分步处理。先把延迟拿掉，画出有理部分的渐近线；再用精确bode图把延迟的影响叠加进去看。因为在幅频特性里，延迟环节的幅值恒为1，只影响相位，所以幅频渐近线不受延迟影响。如果你的重点本来就在幅频上，直接把延迟去掉也无妨。
5.2 转折频率为0导致脚本报错
转折频率为0的情况一般来自原点处的零点或极点。s=0这个因子不会产生转折频率，它只影响初始斜率。我在代码里用abs(real(...)) > eps0的阈值把原点处的零极点单独分离出去，就不会把它们误算成转折频率。
但用户如果手动指定频率范围时把wmin设为0，log10(0)会直接报错。脚本里加了wmin<=0时自动修正到1e-4的兜底逻辑。实际使用中，对数频率轴的起点永远不能是0，最小也得是个很小的正数。
5.3 复根斜率翻倍或曲线断裂
复根斜率翻倍是最隐蔽的bug。如果你自己改过脚本，把共轭复根当成两个独立根处理，那么一个复极点对在转折频率处会贡献-80dB/dec而不是-40dB/dec。曲线断裂则是另一种情况，通常发生在手动删除了wc数组中某个元素但delta没同步删除时。
我用取上半平面的方式从根源上避免翻倍问题。至于断裂，我的建议是不要手动维护wc和delta这两个数组，所有修改都走脚本里的收集、合并流程，除非你非常清楚自己在干什么。
5.4 欠阻尼峰值修正到底怎么做
脚本画的是渐近线，所以在欠阻尼复极点附近，它与实际幅频曲线会有一个偏差。修正方法是：找到对应转折频率wn，算出阻尼比ζ，实际曲线在wn处比渐近线高出20lg(2ζ)分贝。ζ越小，峰值越高。
例如系统G=4/(s^2+s+4)，wn=2，ζ=0.25，修正量20lg(0.5)约-6dB，注意这里公式是20lg(2ζ)，2ζ=0.5，取绝对值得6dB峰值。也就是说，实际曲线在wn附近最高点比渐近线高6dB。补这个峰值时，从渐近线的折点向上画一个圆滑的凸起即可。如果是共轭零点，方向相反，往下凹。
5.5 频率范围选择技巧与图形异常
自动范围一般够用，但遇到转折频率分布跨了好几个数量级的情况，比如一个转折在0.01，另一个在1000，自动范围会变成0.001到10000，中间的线段可能因为图幅限制看不清。这时候手动指定范围更稳妥，优先保证你关心的频段占比大。
图形异常还有一种常见来源：系统增益K为负值。幅频特性取的是|K|，脚本里用abs(k)处理了，所以负增益也能画出正确的幅频渐近线。相位部分就另说了，脚本不画相位，这一点要记住。
最后再分享一个小技巧。想验证脚本画的渐近线对不对，可以把精确幅频曲线叠上去对比：
[mag, ~, wout] = bode(G, freqPts); hold on; semilogx(wout, 20*log10(squeeze(mag)), 'r--');精确曲线应该贴着渐近线走，在转折频率附近圆滑过渡，高低频段完全重合。如果偏差大，基本就是零极点提取环节出问题了。这个交叉验证方法我每次改完脚本都会跑一遍，强烈建议你也试一次。
