简介：面向无线通信、信号处理方向的学习者与研究人员，这份资源聚焦多输入多输出（MIMO）系统的广义奈奎斯特曲线绘制，用MATLAB脚本演示如何由开环传递函数矩阵计算稳定边界，并在复平面上呈现多输入多输出之间的频率响应与相位关系，帮助读者判断系统稳定性、理解信道容量分析中的关键判据。压缩包内共7个文件，均为m脚本，涵盖MIMO广义奈奎斯特绘制、符号传递函数与状态空间转换、模型降阶及示例调用等模块，整体约5KB，体积轻量但环节完整，便于逐步阅读源码并结合理论对照曲线含义。资源已有3166人学习，说明其在课程实验与课题验证中有一定参考价值。对于需要快速搭建MIMO稳定性分析流程、复用绘图脚本或补充模型转换工具的读者，可直接在MATLAB中运行示例，修改系统参数或模拟不同信道环境，观察曲线变化，从而加深对广义奈奎斯特判据与MIMO系统性能特点的理解。

1. 从一次“单回路都稳、整机却自激”说起：MIMO 广义奈奎斯特曲线到底画的是什么

调试多输入多输出（MIMO）回路时，最容易出错的一幕是：把每个回路单独做阶跃，超调、上升时间、相位裕度都漂亮，一合上交叉通道，输出就开始等幅振荡。原因不复杂——单回路裕度只看对角元，耦合项被当成了不存在。要判断这类系统的闭环稳定性，得把整个回路矩阵 L(s)=G(s)K(s) 当作一个整体，看它的特征值 λi(jω) 在复平面上扫出的轨迹离 -1 点多远、净绕了几圈，这就是广义奈奎斯特曲线，也叫特征轨迹。用 MATLAB 把这条曲线画出来本身不难，难的是画对、把绕数数对、让曲线不断线，这几处基本都得自己写。

2. 判据与技术选型：特征轨迹为什么替代不了逐回路裕度

2.1 广义奈奎斯特判据到底在判什么

单回路那套 Nyquist 判据的推导起点是开环传递函数 1+L(s)，多变量系统里对应的量是 det(I+L(s))。把行列式展开成特征值形式，det(I+L) = ∏(1+λi)，其中 λi 是 L(jω) 的特征值。于是闭环特征方程的零点、开环 RHP 极点、以及 det(I+L(jω)) 对原点的绕数之间，仍然满足辐角原理那套关系，只不过“对 -1 点的绕数”要按每条特征轨迹分别统计再求和：闭环稳定 ⟺ Σ Ni = P，其中 Ni 是第 i 条特征轨迹对 -1 点的净绕数（逆时针为正），P 是 L(s) 在右半平面的极点数（重极点按重数计）。

这条判据有两个容易被忽略的前提。一是 L(s) 必须是方阵，即输入数等于输出数；非方阵系统得先做 squared-down 或者直接对 det(I+L) 求绕数，但后者只能给出稳定与否，给不出“离不稳定还有多远”的几何直觉。二是 P 要按 L(s) 的极点算，如果 G 和 K 之间存在零极点对消，而 det(I+L) 里这部分被约掉了，P 的取值就会算错，判据直接失效——我见过一次，模型是从辨识结果里直接拼的，分子分母有几乎相同的因子，P 算出 2，实际闭环是稳的，图形上也看不出异常，排查了大半天。

真正让工程上愿意用特征轨迹的原因，是它把“稳定裕度”变成了看得见的几何量：每条曲线离 -1 点的最近距离 dmin、曲线穿过单位圆时的频率、曲线在 -1 点附近的走向，都能直接读出来。当某条轨迹贴着 -1 点掠过时，你就知道这个方向上的增益再抬一点就会出事，而不像逐回路裕度那样，四个回路都显示 6 dB 裕度、合起来却只剩 1 dB。

2.2 MATLAB 里三条可选路线与取舍

在 MATLAB 里做多变量频域分析，能走的路不止一条，差别在于你要的是图、是数、还是裕度指标。

第一条路线是这次的重点，因为只有它能拿到每条轨迹的原始坐标，后面做配对、着色、自动判据全都依赖这些中间量。第二条路线的价值在于交叉验证：det(I+L) 的绕数应当等于各特征轨迹绕数之和，两个数对不上，基本可以断定配对逻辑或者频率网格出了问题。第三条路线属于“结论层”，diskmargin会返回一个盘型裕度，它考虑的是所有回路同时、同方向地摄动增益和相位，这个指标比逐回路裕度保守也更贴近实际，但它是黑匣子式的输出，出问题时没法告诉你问题出在哪个频率、哪个方向。

这里必须点破一个高频误解：直接对 MIMO 系统调用nyquist(L)，MATLAB 画出来的是每个输入输出对的元素级 Nyquist 曲线，ny×nu 张图，不是特征轨迹。这两者名字像、长得也像，但判据完全不同——元素级曲线对 -1 的绕数在多变量系统里没有稳定性含义。我带的实习生第一次做这个，就是拿元素级图去解释自激现象，结论全反了。

3. 工程骨架：一套可复用的绘图工具的目录与数据流

3.1 目录结构

一次性的脚本能出图，但改两次参数之后就没人敢动了。我习惯把它拆成一个带包目录的小工程，函数名统一挂在+mnyq命名空间下，调用时写mnyq.charLoci(...)，避免和控制工具箱里的函数撞名。

mimo_nyq/ ├── main_demo.m # 一键跑通：建模 -> 算轨迹 -> 画图 -> 打印绕数 ├── +mnyq/ # 包目录，所有算法函数挂在这个命名空间下 │ ├── charLoci.m # L(jw) 逐点特征分解，返回特征值与特征向量 │ ├── pairLoci.m # 特征值配对，保证曲线跨频点连续 │ ├── windingNumber.m # 绕数统计 + 到 -1 点的最小距离 │ ├── plotLoci.m # 频率着色绘图、临界点与单位圆标注 │ └── gainScan.m # 增益扫描，输出稳定域表格 ├── models/ │ ├── model2x2.m # 2x2 强耦合对象，用来验证算法正确性 │ └── model4x4.m # 4x4 过程对象，用来看高维下的"糊图"问题 ├── config/ │ └── plotCfg.m # 频率范围、点数、配色、坐标范围等集中配置 └── tests/ └── test_pairing.m # 人造交叉特征值，专门验证配对逻辑

把模型单独抽成函数，是因为判据验证需要反复换对象；把配置集中到一个 cfg 结构体，是因为频率范围和点数这两个参数几乎每次调试都要改，散落在脚本里迟早会改漏一处。tests/test_pairing.m看起来有点小题大做，但配对逻辑是全工程最容易出隐性错误的地方——它错了图还是能画出来，只是曲线偷偷换了道，肉眼很难发现。

3.2 从模型到曲线的六步数据流

整条链路的中间变量形状变化很容易记混，尤其是freqresp返回的三维数组，索引顺序踩过一次坑就会记一辈子。

数据流里唯一有分支的是最后一步：要图就走plotLoci，要结论就走windingNumber，两个函数吃的是同一份配对后的lam。这种设计的好处是，图上的几何直觉和判据的数值结论永远来自同一批数据，不会出现“图看着稳、代码判不稳”的对不上账。

4. 核心代码拆解：从 L(jω) 到能看懂的特征轨迹

4.1 建模与频率网格：先把网格铺对

% 2x2 被控对象 + 对角控制器，故意保留较强的交叉通道 s = tf('s'); G = [ 1/(s+1), 0.6/(s+2); 0.4/(s+3), 2/(s+0.5) ]; K = diag([1.5, 0.8]); % 先只看比例控制，耦合影响最明显 L = G * K; % 回路传递矩阵，注意是矩阵乘法 L = minreal(L, 1e-6); % 消去数值上的零极点对，避免响应毛刺 % 频率网格：对数均匀铺底 + 在 0.1~10 rad/s 区间加密 w = [ logspace(-3, -1, 80), linspace(0.1, 10, 400), logspace(1, 4, 120) ]; w = unique(w(:)); % 去重并转成列向量

minreal那一步不是可选项。辨识或者手工拼接出来的模型，分子分母常带 1e-8 量级的近似公因子，freqresp在少数频点上会算出明显偏离的值，反映到特征轨迹上就是曲线上突然冒出一个尖刺，看图的人会以为系统在某个频率有谐振。容差 1e-6 是经验值，太小了消不掉，太大了会把真实的相近极点也吃掉。

频率网格用了三段拼接而不是单个logspace。原因是特征轨迹穿越 -1 点附近的区域通常集中在穿越频率周围，对数均匀采样在这段区间点数偏少，绕数统计容易漏掉一次小回环。低频段用 80 点是因为低频特征值变化平缓，不需要密；高频段用logspace(1,4,120)是因为曲线在高频会收敛到常值矩阵的特征值，点太密只是浪费计算。真要省事，logspace(-3,4,2000)也能用，代价是eig被调用两千次，4x4 以上的系统会明显卡顿。

4.2 逐点特征分解：freqresp 返回的三维数组怎么用

function [lam, V, w] = charLoci(L, w) %CHARLOCI 计算 L(jw) 的特征值（特征轨迹）与对应特征向量 % lam : nw x n 复数矩阵，第 k 行是频率 w(k) 处的 n 个特征值 % V : n x n x nw，第 k 页的列向量是 lam(k,:) 对应的特征向量 n = size(L, 1); nw = numel(w); fr = freqresp(L, w); % 返回 n x n x nw 的三维数组 lam = zeros(nw, n); V = zeros(n, n, nw); for k = 1:nw [Vk, D] = eig(fr(:,:,k)); % 单频点特征分解 lam(k,:) = diag(D).'; % 转置成行，方便按频率堆叠 V(:,:,k) = Vk; end end

freqresp的返回顺序是“输出×输入×频率”，第三个维度才是频率，这一点和bode返回的mag维度一致但和很多人的直觉相反。写成fr(:,:,k)取第 k 个频点的 n×n 复矩阵，直接丢给eig即可。

必须同时保留特征向量V，这是后面配对的前提。eig返回的特征值顺序是算法内部决定的，同一个系统在相邻两个频点上，第 1 个特征值可能属于完全不同的分支——eig本身不保证任何连续性。如果只留特征值，画出来的曲线会在某些频点突然“换道”，看起来像是断成了几截，实际上数据没错，是索引错了。

顺带说一个数值细节：eig对非对称实矩阵返回的是复数特征值，理论上 L(jω) 在实频率上的特征值不保证共轭成对出现（除非 L 是实系数矩阵在实频率取值，此时才共轭成对）。所以判断“某条轨迹是否穿到右半平面”不能简单看虚部符号，得看实部。

4.3 特征值配对：让曲线不断线的关键一步

function [lam, V] = pairLoci(lam, V) %PAIRLOCI 按特征向量相关性把相邻频点的特征值接上，消除顺序跳变 % 思路：v_i(jw_k) 与 v_j(jw_{k-1}) 的内积模越大，越可能是同一条分支 [nw, n] = size(lam); for k = 2:nw C = abs(V(:,:,k)' * V(:,:,k-1)); % n x n，C(i,j)=|<v_i^k, v_j^{k-1}>| usedR = false(n,1); usedC = false(n,1); perm = 1:n; for t = 1:n Cc = C; Cc(usedR,:) = -inf; Cc(:,usedC) = -inf; % 已占用的行列屏蔽 [~, idx] = max(Cc(:)); [r, c] = ind2sub(size(Cc), idx); % r: 新频点序号, c: 旧频点槽位 perm(c) = r; % 槽位 c 放第 r 个特征值 usedR(r) = true; usedC(c) = true; end lam(k,:) = lam(k, perm); V(:,:,k) = V(:, perm, k); end end

这是整套代码里最值得单独写一个函数的地方。核心假设是：频率变化足够小时，同一条特征轨迹的特征向量方向也只会小幅旋转，所以用相邻频点特征向量的内积模作为“相似度”，做一次贪心指派就能把顺序对齐。用内积模而不是实部，是因为特征向量带任意相位，v和e^{jθ}v表示同一个方向，取模消掉这个自由度。

Cc(usedR,:) = -inf这几行是在做互斥：一旦某个新序号被分配出去，它所在的整行和对应槽位的整列都要屏蔽，否则贪心会重复指派。贪心不保证全局最优，遇到两条轨迹在某个频点几乎重合（特征值简并）时可能选错，但代价是错一次只影响局部一段，不会像排序法那样整条曲线重排。要更稳可以换成匈牙利算法matchpairs(-C, 1e3)，R2019a 之后才有，内部实现是全局最优指派，代价是每步 O(n³)，n 大于 8 之后比贪心慢一个数量级。

配对做完之后还有个隐性收益：V的顺序也一起对齐了，后面要算特征向量的相位变化率、或者做模态分析，直接就能用。

4.4 绕数与临界点：把“看图”变成“判据”

function [N, dmin] = windingNumber(lam, w, target) %WINDINGNUMBER 统计每条特征轨迹绕 target（默认 -1）的净圈数 % N : 1 x n，逆时针为正、顺时针为负的净绕数 % dmin : 所有轨迹到 target 的最小距离，衡量稳定裕度 if nargin < 3, target = -1; end [nw, n] = size(lam); N = zeros(1, n); dmin = inf; for i = 1:n z = lam(:,i) - target; % 平移到以临界点为原点 d = abs(z); dmin = min(dmin, min(d)); dth = diff(unwrap(angle(z))); % 解卷绕后的相位增量 N(i) = round(sum(dth) / (2*pi)); % 净绕数 = 总相位增量 / 2pi end end

unwrap是这里的命门。直接对angle的结果做差分，角度在 ±π 处会跳变 2π，累加出来的总相位增量就完全不对了。unwrap把跳变补回去，sum(dth)/(2π)才是真实的净绕数；round用来吸收浮点累积误差，正常情况下结果应该非常接近整数，如果算出来是 0.5 这种值，说明频率网格太粗，轨迹在某个频段绕了半圈就被跳过去了。

还有一个绕数统计的范围问题：上面这段只处理了 ω 从 0 到 +∞ 的半条轨迹。完整的 Nyquist 围线包含负频率部分，而 λ(-jω) 与 λ(jω) 共轭，所以完整绕数等于半条轨迹绕数的两倍，再加上无穷远处的圆弧贡献。对严格真有理的 L(s)（分母阶次高于分子），ω→∞ 时 L→0，圆弧收缩成一个点，不贡献绕数；对双真系统，L(∞) 是个常值矩阵，那一段得按常值矩阵的特征值单独补一段路径。代码里没写这段是因为示例对象都是严格真的，换对象时必须自己确认这一点，不然判据会差一圈。

dmin这个输出用起来比绕数更顺手。绕数只有 0、±1、±2 几个离散值，而dmin是连续的，比值小于 0.3 就说明曲线贴着临界点走，哪怕当前还是稳定的，增益稍微一动就会翻车。

4.5 绘图：频率着色与坐标等比例

function plotLoci(lam, w, cfg) %PLOTLOCI 以频率着色的方式画特征轨迹，并标注临界点与单位圆 figure('Color','w'); hold on; cmap = parula(numel(w)); % 新版可换 turbo，对比度更高 for i = 1:size(lam,2) for k = 1:numel(w)-1 % 逐段着色，才能体现频率方向 plot(real(lam(k:k+1,i)), imag(lam(k:k+1,i)), '-', ... 'Color', cmap(k,:), 'LineWidth', 1.4); end plot(real(lam(1,i)), imag(lam(1,i)), 'o', 'MarkerSize',5); % 起点 plot(real(lam(end,i)), imag(lam(end,i)), 's', 'MarkerSize',5); % 高频端 end plot(-1, 0, 'k+', 'MarkerSize', 14, 'LineWidth', 2); % 临界点 th = linspace(0, 2*pi, 400); plot(cos(th), sin(th), 'k:'); % 单位圆参考 axis equal; grid on; xlabel('Re'); ylabel('Im'); caxis(log10([w(1) w(end)])); colorbar; title(sprintf('广义奈奎斯特曲线：%d 条特征轨迹', size(lam,2))); end

axis equal是必须的。默认坐标轴会按数据范围自动拉伸，一条本来离 -1 点还有 0.5 距离的曲线，在横向被压扁之后看起来像是直接穿过去了，看图的人会得出完全相反的结论。加上axis equal，实轴虚轴同比例，几何距离才有意义。

逐段plot而不是一次性画整条，是为了让颜色随频率变化——对数频率下曲线的疏密本身就携带信息，低频段点稀、高频段点密，颜色一叠上去，哪一段是穿越区一眼就能认出来。缺点是段数多的时候渲染慢，2000 个频点、4 条轨迹就是 8000 次plot调用，出图能等上好几秒；要批量出图就换成surface或者patch一次性铺，牺牲一点可控性换速度。

caxis在较新的版本里改名叫clim，老代码搬过去会报错，改成clim(log10([w(1) w(end)]))即可。

5. 避坑与排查：五条踩过的坑

5.1 现象：曲线莫名其妙断成几截，或者两条轨迹在某个频率交换了位置

原因是eig不保证相邻频点的特征值顺序一致，而排序法（按实部或按相角排序）在特征值接近简并时会判错，两条轨迹整段换道，图上看起来就像断了。

解决办法就是 4.3 节的向量配对，用特征向量相似度做指派。判断有没有配对成功，有个简单办法：把lam每列相邻点的差值取模，如果某处出现远大于平均值的跳变（比如比中位数大两个数量级），说明那里没接上。tests/test_pairing.m就是干这个的，构造一对交叉的特征值函数，跑一遍看配对结果和解析解是否一致。

5.2 现象：绕数算出来不是整数，或者和 det(I+L) 的绕数对不上

两个原因。一是频率网格不够密，轨迹在穿越区绕的小回环被采样的直线段“切掉”了；二是漏算了负频率和无穷远弧段的贡献。

处理办法分两步：先用det(eye(n)+L(jw))单独算一遍总绕数作为真值，再拿特征轨迹的绕数之和去比；对不上就加密网格，特别是穿越频率附近，用linspace局部补点。我在 4x4 对象上遇到过差一圈的情况，最后发现是高频段logspace(1,4,120)太稀，一条轨迹在 500 rad/s 附近有个很小的回环被完全跳过，补到 400 点之后才对上。

5.3 现象：图上四条轨迹挤成一团，完全看不出哪条离 -1 最近

高维系统常见。n 大于 3 之后，特征轨迹在复平面上的分布会非常拥挤，尤其是对象各通道增益差异大的时候，有的轨迹贴近原点，有的跑到几十以外，坐标范围一拉开，临界点附近全糊了。

解决办法是先做通道缩放，把 L 的输入输出按稳态增益归一化，让各通道量级接近，再画图；实在挤不开就分面画，subplot每张图只放一两条轨迹，配上各自的dmin标注。缩放的意义不只是好看——广义奈奎斯特判据本身对相似变换不变，缩放不改变绕数，所以可以放心做。

5.4 现象：把nyquist(L)的图当成广义奈奎斯特曲线

这是最典型的一个。nyquist(L)对 MIMO 系统画的是元素级曲线，2x2 系统出四张子图，每张对应一个输入输出对，判据和特征轨迹完全不是一回事。

区分方法很直接：元素级图的曲线数量是 ny×nu，特征轨迹的数量恒等于 n；再不行就看代码里有没有eig(freqresp(...))这一步，没有就不是。写报告的时候我一般会两张图都放，元素级用来说明耦合强度，特征轨迹用来下稳定性结论，读者不会混淆。

5.5 现象：模型改了一个参数，判据从稳定变成不稳定，但阶跃响应看着正常

多半是零极点对消导致 P 算错。L = G*K之后如果没做minreal，分子分母里残留的近似公因子会让pole(L)返回一堆虚假的右半平面极点，P 从 0 变成 2，判据立刻翻转。

排查顺序：先zpk(L)看一眼零极点分布，有没有相距小于 1e-4 的对子；再做minreal(L, 1e-6)；最后用闭环极点max(real(pole(feedback(L, eye(n)))))交叉验证。闭环极点全在左半平面、而判据说不行，那一定是 P 算错了，不要反过来怀疑判据。

提示：判据代码和闭环极点验证最好永远成对出现。绕数只是个数字，闭环极点才是最终事实，两者对不上时优先信极点，然后回头查频率网格、配对和 P 的计算。

6. 二次开发：从一张图到一个能批量判稳的小工具

单张图解决的是“看一眼”的需求。实际项目里更常见的场景是：控制器结构定了，要扫一遍增益看稳定边界在哪；或者对象有多个工况，每个工况都要判一次。这时候把 4.4 节的绕数函数往外扩一层就行。

function T = gainScan(G, K0, gains, w) %GAINSCAN 扫描标量增益 k，用广义奈奎斯特判据判断闭环稳定性 n = size(G,1); T = table(); for k = gains L = minreal(G * (k*K0), 1e-6); [lam, V] = mnyq.charLoci(L, w); [lam, ~] = mnyq.pairLoci(lam, V); N = mnyq.windingNumber(lam, w, -1); P = sum(real(pole(L)) > 0); % 开环 RHP 极点数 Ntot = sum(N); % 各轨迹绕数之和 Tcl = feedback(L, eye(n)); poleOK = all(real(pole(Tcl)) < -1e-8); % 闭环极点真值 T = [T; table(k, P, Ntot, Ntot == P, poleOK)]; end end

关键在最后两列：Ntot == P是判据给出的结论，poleOK是闭环极点的真实结论，两列并排放在表里，跑十次只要有一次不一致，立刻就能定位到是哪个增益点出了问题。我一般还会在表里加一列min(abs(lam(:)+1))，记录该增益下离临界点的最近距离，稳定域边界附近这一列会趋近于零，比单纯看 true/false 有信息量得多。

进一步可以做的扩展有三个方向。一是把w做成自适应网格：先粗扫一遍，找出dmin的极小值所在频段，再在那个频段局部加密重算，能省掉一大半计算量。二是和diskmargin交叉验证，两者的结论趋势应当一致，diskmargin偏保守，如果出现盘裕度说稳定而绕数说不稳定，那基本可以确定是 P 或网格的问题。三是出报告：把特征轨迹图、绕数表、稳定边界图拼成一张，用exportgraphics直接导出 PDF，评审的时候比一堆散图省事。

有个细节值得单独说：gainScan里每次循环都重建了L，对tf对象来说这步很慢。如果G和K0都是状态空间模型，把G和K0提前转成ss、用k直接缩放K0的 B、C 矩阵，循环速度能快一个数量级；4x4 对象扫 50 个增益点，tf版本要跑两分多钟，状态空间版本十几秒就完事。

踩过最疼的一次是：扫增益时忘了在循环里做minreal，某个增益点恰好让两个极点靠得极近，pole判定翻车，判据给出的稳定域边界偏了整整 0.4。从那以后我每次写这类批量扫描，都强制在循环体内做一次模型化简，并且在结果表里保留闭环极点验证那一列——宁可多算一遍，也不想再拿着一张边界位置错误的图去开会。希望这些能帮到你，少走一遍我走过的弯路。

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