简介：面向无线通信、信号处理方向的学习者与研究人员，这份MATLAB代码资源聚焦多输入多输出（MIMO）系统的广义奈奎斯特曲线绘制，帮助评估系统稳定性与频率响应特性。压缩包共7个文件，全部为.m脚本，整体约5KB，涵盖传递函数符号化转换、模型降阶、广义奈奎斯特绘图及示例主程序等模块，便于按流程理解从建模到成图的完整链路。已有3163人学习下载，说明其在教学与工程分析中具备一定参考价值。读者可通过运行示例脚本，观察复数传递函数矩阵在复平面上的稳定边界表现，并结合参数调整与信道环境模拟，加深对MIMO系统相位特性与稳定条件的理解，同时提升MATLAB编程与通信理论结合应用的能力。

1. 从一张“绕圈”的曲线说起：MIMO 广义奈奎斯特到底在画什么

做多变量反馈设计的同行大概都遇到过这种场面：单回路里看奈奎斯特曲线，一眼就能判断稳定裕度，可一旦系统变成两输入两输出，把每个通道的奈奎斯特图分别画出来，曲线看着都离临界点很远，闭环却偏偏发散。问题不在图画错了，而在于多变量系统的稳定性不能靠“逐个通道看”来判断——回路之间存在耦合，一个通道的增益变化会牵动另一个通道的等效开环特性。广义奈奎斯特曲线就是为解决这件事而生的：它把 MIMO 开环传递矩阵 L(s) 在每个频率点的特征值 λ_i(jω) 取出来，把这些特征值轨迹画在同一张复平面上，用它们绕 (-1, 0) 点的圈数来判断闭环稳定性。这份 MATLAB 资源围绕的正是这条曲线的绘制流程，适合正在做多变量频域设计、需要把稳定性判据落到图上的控制方向开发者。

2. 技术底座与目录结构：这套 MATLAB 工程由哪些文件撑起来

2.1 为什么选频域特征值轨迹而不是时域仿真

时域仿真能告诉你“这一组参数下发散”，但说不清“离发散还有多远、往哪个方向调参数更安全”。广义奈奎斯特判据的价值在于给出裕度信息：特征值轨迹离临界点的最近距离、绕圈方向、穿越负实轴的频率区间，这些都是可以直接指导参数整定的量。MATLAB 在这个场景下的优势是eig、freqresp、nyquistplot这些函数现成，传递矩阵可以直接用tf或ss对象拼出来，不用自己写频响计算。常见做法是把被控对象写成状态空间形式，控制器写成传递矩阵，两者相乘得到开环 L(s)，再在关心的频段上扫点求特征值。

需要提前说清一个概念边界：广义奈奎斯特判据针对的是方阵开环传递矩阵，输入输出维数必须相等。如果你的系统是 2 入 3 出这类非方情形，得先做输入或输出变换把它方阵化，否则特征值无从谈起。这一点在动手写代码前就要确认，不然后面所有图都是错的。

2.2 工程目录的典型组织方式

这类绘制工程一般不会只有一个脚本，按职责拆开更利于二次开发。下面是我拆这类项目时常见的目录骨架，文件名按功能命名，你可以对照自己手上的资源包核对：

mimo_nyquist/ ├── main_nyquist.m % 主入口：定义对象、控制器，调用绘制流程 ├── model/ │ ├── build_plant.m % 构造被控对象传递矩阵 G(s) │ └── build_controller.m % 构造控制器矩阵 K(s) ├── core/ │ ├── open_loop.m % 计算开环 L(s)=G(s)K(s) │ ├── eig_sweep.m % 频段扫描 + 逐点特征值分解 │ └── winding_number.m % 绕圈数统计与临界点距离 ├── plot/ │ ├── plot_gnq.m % 广义奈奎斯特曲线绘制 │ └── plot_margin.m % 裕度标注（最近距离、穿越频率） └── config/ └── freq_range.m % 频率范围与采样密度配置

这个划分的核心思路是“模型构造 / 数值计算 / 绘图”三层分离。模型层只负责把 G(s)、K(s) 拼出来，计算层只吃传递矩阵吐特征值序列，绘图层只吃数据画图。好处是换一个被控对象时只改model/下的两个文件，绘图和绕圈数统计完全复用；反过来想换绘图风格，也不碰数值代码。我一般会把频率范围和采样点数单独抽到config/，因为这两个参数对曲线形状影响极大，集中管理比散落在脚本里好排查。

2.3 关键函数依赖与版本注意

整个流程依赖 Control System Toolbox，freqresp、tf、ss、eig都在这个工具箱里。如果你的 MATLAB 没装这个工具箱，tf和freqresp会直接报未定义函数，这是新手最容易卡住的第一步。另外nyquistplot在较新版本里对 MIMO 对象的支持有变化，自己用plot画特征值轨迹反而更可控——这也是为什么这类工程通常不直接调nyquistplot，而是手动算特征值再plot。

3. 核心计算链路：从开环矩阵到特征值轨迹的逐段实现

3.1 构造开环传递矩阵 L(s)

一切从开环开始。设被控对象 G(s) 是 2×2，控制器 K(s) 也是 2×2，开环 L(s) = G(s)K(s)。用 MATLAB 的传递矩阵写法，每个元素是一个tf对象，拼成矩阵后矩阵乘法自动按传递函数运算：

% build_plant.m —— 构造 2x2 被控对象传递矩阵 s = tf('s'); % 对角通道：两个一阶惯性环节，时间常数不同 g11 = 1 / (s + 1); g22 = 2 / (0.5*s + 1); % 耦合通道：交叉项，体现回路间相互影响 g12 = 0.3 / (s + 2); g21 = 0.2 / (s + 3); G = [g11, g12; g21, g22]; % 控制器：简单的比例 + 积分，对角结构 k11 = 1.5 + 0.8/s; k22 = 1.2 + 0.5/s; K = [k11, 0; 0, k22]; % 开环传递矩阵 L = G * K;

这段代码的关键点有三个。第一，s = tf('s')声明了传递函数变量，后面所有1/(s+1)这类写法才会被识别为传递函数而不是数值除法，漏了这一行会得到一堆标量而不是tf对象。第二，交叉项g12、g21不能省，它们正是“单通道看着稳、整体却发散”的根源，去掉耦合项再画广义奈奎斯特，你会得到一个过于乐观的结论。第三，控制器用对角结构是为了先验证基础流程，实际工程里 K 可能是满阵，矩阵乘法G*K照样成立，不需要改代码结构。

参数上，时间常数和增益决定了特征值轨迹的走向。g11的极点 -1 和g22的极点 -2 会让低频段特征值偏向实轴正方向，耦合项系数 0.3、0.2 越大，两条特征值轨迹分得越开。调这几个数，你能直观看到曲线形状怎么变，这是理解判据最快的办法。

3.2 频段扫描与逐点特征值分解

有了 L(s)，下一步是在对数分布的频率点上求每个点的特征值。频率范围要覆盖系统的主要动态：下限低于最慢极点一个数量级，上限高于最快极点一个数量级，否则曲线两端会“缺头少尾”，绕圈数统计算不准。

% eig_sweep.m —— 频段扫描并逐点求特征值 function [lambda_all, w] = eig_sweep(L, wmin, wmax, npts) % 对数分布频率点，避免低频段采样过疏 w = logspace(log10(wmin), log10(wmax), npts); % 计算开环频率响应，返回 2x2xnpts 的三维数组 Lresp = freqresp(L, w); nout = size(Lresp, 1); lambda_all = zeros(nout, numel(w)); % 每列存一个频率点的全部特征值 for k = 1:numel(w) Lk = Lresp(:, :, k); % 取出第 k 个频率点的复矩阵 lam = eig(Lk); % 特征值分解 % 按实部排序，保证相邻频率点的特征值轨迹连续 [~, idx] = sort(real(lam)); lambda_all(:, k) = lam(idx); end end

这里有个血泪经验：eig返回的特征值顺序是任意的，不同频率点之间顺序可能跳变，直接画出来会得到一堆交叉乱线，根本看不出轨迹。解决办法是按实部（或按与上一频率点特征值的最近距离）排序，让同一支轨迹在相邻频率点上“接得上”。上面用实部排序是最简单的做法，对大多数工程对象够用；如果两支轨迹实部很接近、排序仍会跳，就改成“每个特征值匹配上一列中距离最近的那个”，代码多几行但更稳。

freqresp返回的是三维数组，第一维输出、第二维输入、第三维频率点，取Lresp(:,:,k)才是第 k 个频率点的复矩阵。这个索引顺序搞反是常见翻车点，症状是eig报维度不匹配或者算出来的特征值完全不对。npts一般取 500 到 2000，点太少曲线折线感重、绕圈数统计会漏，点太多绘图卡顿但精度提升有限，1000 是个稳妥的起点。

3.3 绕圈数统计与临界点距离

广义奈奎斯特判据的核心是数特征值轨迹绕 (-1, 0) 的净圈数。数值上不能靠“看”，得用辐角累积法：把每个频率点相对临界点的复向量角度算出来，相邻点角度差做解缠绕，累加后除以 2π 就是绕圈数。

% winding_number.m —— 统计单支特征值轨迹绕 (-1,0) 的净圈数 function N = winding_number(lam_track, crit) if nargin < 2 crit = -1 + 0i; % 默认临界点 -1+j0 end % 相对临界点的复向量 rel = lam_track - crit; % 相邻点角度差，unwrap 处理 ±pi 跳变 dtheta = diff(unwrap(angle(rel))); % 净绕圈数 = 总角度变化 / 2π N = sum(dtheta) / (2*pi); end

unwrap这一步不能省。angle返回的是 (-π, π] 区间内的值，轨迹穿过负实轴时角度会从 π 跳到 -π，不 unwrap 就会把一次正常穿越算成反向绕圈，圈数直接错。crit参数留成可配置，是因为有些教材把临界点定义在 (-1, 0)，有些在原点附近做变换，留个口子方便对照。

除了圈数，实际调参更关心“离临界点多近”。在eig_sweep的结果上再算一个最小距离：

% 所有特征值轨迹到临界点的最小距离，作为稳定裕度的粗略指标 all_lam = lambda_all(:); min_dist = min(abs(all_lam - (-1 + 0i)));

这个距离不是严格的稳定裕度定义（多变量系统有更严谨的 μ 分析），但作为快速筛查足够用：距离越小，参数越接近危险区。我一般会在绘图时把这个最近点和对应频率标出来，调参时盯着它走。

4. 绘图与结果解读：把特征值轨迹画成能读的图

4.1 广义奈奎斯特曲线绘制

数据齐了，绘图本身不复杂，难点在“让图能读”：多支轨迹要区分颜色、临界点要标出、单位圆和负实轴要画参考线。

% plot_gnq.m —— 绘制广义奈奎斯特曲线 function plot_gnq(lambda_all, w) figure; hold on; grid on; ntrk = size(lambda_all, 1); colors = lines(ntrk); % 每支轨迹一种颜色 for i = 1:ntrk plot(real(lambda_all(i,:)), imag(lambda_all(i,:)), ... 'Color', colors(i,:), 'LineWidth', 1.2, ... 'DisplayName', sprintf('\\lambda_%d', i)); % 标出低频起点，方便判断轨迹走向 plot(real(lambda_all(i,1)), imag(lambda_all(i,1)), ... 'o', 'Color', colors(i,:), 'MarkerFaceColor', colors(i,:)); end % 临界点 plot(-1, 0, 'rx', 'MarkerSize', 12, 'LineWidth', 2, ... 'DisplayName', 'critical point (-1,0)'); % 单位圆参考 theta = linspace(0, 2*pi, 200); plot(cos(theta), sin(theta), 'k--', 'LineWidth', 0.5); axis equal; xlabel('Real'); ylabel('Imag'); title('Generalized Nyquist Plot (eigenvalue loci)'); legend('Location', 'best'); end

axis equal必须加，否则实轴虚轴比例不一致，圆看起来是椭圆，绕圈关系会被视觉扭曲。低频起点用实心圆点标出来，是因为轨迹方向（随频率增大往哪走）对判断绕圈方向很关键，光看一条线分不清头尾。lines(ntrk)生成区分度足够的颜色序列，轨迹支数多时比手动指定颜色省事。

4.2 结果怎么读：三个判断动作

图出来之后，按顺序做三件事。第一，看每支轨迹是否包围 (-1, 0)：包围则对应闭环不稳定（在开环无右半平面极点的前提下），这是判据的直接结论。第二，看轨迹离临界点的最近距离和对应频率，距离小说明该频段裕度不足，调参优先动这个频段的增益。第三，看轨迹穿越负实轴的位置，穿越点越靠左（越接近 -1），相位裕度越紧张。

需要提醒的是，广义奈奎斯特判据的完整形式还要考虑开环右半平面极点数 P，闭环不稳定极点数 Z = P + N（N 为绕圈数，符号约定因教材而异）。如果你的对象本身开环不稳定（比如倒立摆类），P 不为零，光看“包不包围”会得出错误结论。这类工程里通常会在winding_number旁边再放一个pole(L)统计右半平面极点数，两个数一起看才完整。

4.3 和单回路奈奎斯特图的对照

把同一系统的单通道奈奎斯特图和广义奈奎斯特图并排看，是理解耦合影响最直观的方式。单通道图里每个对角元素的曲线可能都远离 -1，但广义图里特征值轨迹却绕了进去——这个反差正是多变量系统“通道裕度不等于系统裕度”的证据。做二次开发时，我建议保留一个对照绘图函数，把两种图放一个 figure 的两个 subplot 里，汇报和自查都用得上。

5. 避坑与排查：这类绘制工程最容易翻车的五个地方

5.1 曲线画出来是乱麻，轨迹互相穿插

现象：绘图结果是一团交叉的线，看不出任何连续轨迹。原因：eig返回顺序在相邻频率点间跳变，没有做排序或匹配。解决：在eig_sweep里按实部排序，或者用“最近邻匹配上一列特征值”的方式重排；排序后仍乱，检查频率点是否过疏导致轨迹跳变过大，把npts提到 2000 再试。

5.2 绕圈数算出来是小数或明显不对

现象：winding_number返回 0.5、1.7 这类非整数。原因：angle的 ±π 跳变没做unwrap，或者频率范围没覆盖完整轨迹导致首尾不闭合。解决：确认unwrap在diff之前调用；把频率下限再降一个数量级、上限再升一个数量级，让轨迹两端都收敛到稳定值；检查crit是否和判据约定一致。

5.3 非方阵系统直接报错

现象：eig报“矩阵必须为方阵”。原因：开环 L(s) 是 m×n 且 m≠n，特征值无定义。解决：广义奈奎斯特判据只适用于方阵开环，非方系统要先做输入/输出变换（比如取部分输出做反馈、或引入权重矩阵做方阵化），这一步没有通用公式，得按具体系统结构设计，别硬套。

5.4 换了对象后图完全不对，但代码没报错

现象：曲线形状离谱，和预期动态对不上。原因：freqresp返回的三维数组索引顺序搞反，把输入维当成了频率维；或者tf变量s在函数内部重新声明导致作用域混乱。解决：在eig_sweep里打印size(Lresp)确认是[nout, nin, nfreq]；把s = tf('s')统一放在模型构造函数内部，不要跨函数共享。

5.5 低频段曲线发散到无穷远

现象：轨迹在低频端冲向无穷。原因：控制器含积分环节（1/s），低频增益趋于无穷，特征值模值随之爆炸。解决：这是物理事实不是 bug，绘图时把频率下限设在积分起作用之上一点，或者对幅值做截断显示；判据分析时低频段本来就要单独处理，别指望一张图看全。

6. 二次开发：把绘制流程改造成自己的稳定性检查工具

拿到这套代码后，最值得做的改造不是换绘图配色，而是把它变成一个能批量跑参数的检查函数。我的做法是把main_nyquist.m里的流程包成一个函数，输入是对象参数和控制器参数，输出是绕圈数、最小距离、是否稳定三个量，这样就能套for循环扫参数网格，找出稳定域边界。

% check_stability.m —— 批量稳定性检查的封装 function [is_stable, min_dist, N_total] = check_stability(p) % p 是参数结构体：增益、时间常数、耦合系数 G = build_plant(p); K = build_controller(p); L = G * K; [lam, ~] = eig_sweep(L, 1e-2, 1e3, 1000); % 统计所有轨迹的净绕圈数之和 N_total = 0; for i = 1:size(lam, 1) N_total = N_total + winding_number(lam(i,:)); end % 开环右半平面极点数 P = sum(real(pole(L)) > 0); % 闭环稳定判据：Z = P + N_total，Z 为 0 则稳定 is_stable = (P + round(N_total)) == 0; min_dist = min(abs(lam(:) - (-1 + 0i))); end

这个封装的关键是把“看图判断”变成“数值判断”，round(N_total)处理浮点误差，P从pole(L)直接统计。有了它，扫参数就是两层循环的事，把is_stable为真的参数点收集起来，稳定域边界自然浮现。验证方法也简单：拿一组已知稳定的参数跑一遍，is_stable应为真、min_dist应明显大于零；再故意把增益调大到发散，is_stable应翻假。两个方向都对，封装才算可信。

还有一个容易忽略的改造点：频率范围不该写死。不同对象的动态时间尺度差几个数量级，1e-2到1e3只适合秒级对象。我一般从pole(G)和zero(G)的模值自动推频率范围，取最小极点模值的百分之一到最大极点模值的百倍，这样换对象不用手调。这个自动推范围的逻辑放进config/freq_range.m，和手动配置并存，默认走自动、特殊情况手动覆盖。

从那以后我每次把这类频域判据代码交给别人用之前，都强制拿一个已知答案的简单系统跑一遍对照——单输入单输出退化成经典奈奎斯特，结果必须和手算一致，这个回归测试能挡掉大部分索引和符号约定的低级错误。希望帮到你。

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