简介：面向无线通信、信号处理方向的学习者与研究人员，这份MATLAB代码资源聚焦多输入多输出（MIMO）系统的广义奈奎斯特曲线绘制，用于评估系统稳定性与频率响应特性。资源包共7个文件，全部为.m脚本，压缩包约5KB，涵盖传递函数符号化转换、模型降阶、广义奈奎斯特绘图及示例主程序等模块，便于按流程理解从建模到成图的完整链路。已有3163人学习下载，说明其在教学与工程分析中具有一定参考价值。通过阅读与运行这些脚本，读者可掌握nyqmimo类函数对开环传递函数矩阵的处理方式，观察复平面上的稳定边界，并结合参数调整与信道环境模拟，加深对MIMO系统相位特性与稳定条件的理解，同时提升MATLAB编程与通信系统分析能力。
1. 广义奈奎斯特曲线：为什么 MIMO 系统的稳定性不能靠“一条曲线”判断
如果你在 SISO 系统里调过反馈，大概率形成过一种肌肉记忆：开环传函画一条奈奎斯特曲线，看它围没围住 -1 点，围住就闭环不稳定，没围住就稳。这套判断在单输入单输出里几乎成了条件反射。但把同样的思路搬到多输入多输出 MIMO 系统上，会立刻翻车——你面对的不再是一条曲线，而是一族曲线，围住 -1 点这件事本身都失去了定义。广义奈奎斯特曲线要解决的，正是这个“标量判据失效”的问题：它把稳定性判断从“看一个点被围几次”升级为“看一组特征值轨迹绕点几圈”，用矩阵特征值的轨迹替代单条开环频率响应。搞懂它，你才能理解为什么 MIMO 反馈系统的稳定裕度是个矩阵问题，而不是逐回路调参能凑出来的玄学。
2. 从 SISO 到 MIMO：奈奎斯特判据为什么必须“广义化”
2.1 标量判据的成立前提，在 MIMO 里悄悄消失了
SISO 奈奎斯特判据的底层是辐角原理：开环传函 (L(s)) 在右半平面有 (P) 个极点，闭环稳定的充要条件是奈奎斯特曲线逆时针包围 (-1+j0) 的圈数 (N=P)。这个结论成立依赖一个隐含前提——闭环特征方程是 (1+L(s)=0)，一个标量方程，一个标量函数。到了 MIMO，开环传函变成一个 (m\times m) 的有理函数矩阵 (L(s))，闭环特征方程变成 (\det(I+L(s))=0)。注意，稳定性由行列式决定，而不是由某个回路的 (1+L_{ii}(s)) 决定。行列式是各特征值的乘积，这意味着任何一个特征值穿过 -1 点附近，行列式都可能变号。换句话说，MIMO 的稳定性是“所有特征值共同作用”的结果，逐回路单独看奈奎斯特曲线，等于把耦合项全部丢掉，判断自然不可靠。
常见的实现思路是把行列式判据直接搬上频率轴：令 (s=j\omega)，画 (\det(I+L(j\omega))) 随 (\omega) 变化的轨迹，看它是否包围原点。这确实是一种广义化，但它有个致命弱点——行列式是特征值的乘积，多个特征值同时靠近原点时，乘积的轨迹可能互相抵消，把不稳定的迹象“抹平”。所以更稳妥的做法是直接追踪 (L(j\omega)) 的特征值 (\lambda_i(j\omega))，把每个特征值当成一条独立的复平面轨迹，这就是广义奈奎斯特曲线（generalized Nyquist / characteristic loci）的核心思想。
2.2 特征轨迹判据：把“围点”换成“围 -1 的圈数统计”
设 (L(s)) 是 (m\times m) 开环传函矩阵，(\lambda_i(s)) 是它的第 (i) 个特征值。广义奈奎斯特判据可以写成：
[ \sum_{i=1}^{m} N_i = P ]
其中 (N_i) 是第 (i) 条特征轨迹 (\lambda_i(j\omega)) 逆时针包围 (-1+j0) 的净圈数，(P) 是 (L(s)) 在右半平面的极点数（计重数）。这个式子和 SISO 的 (N=P) 形式几乎一样，但含义变了：不是一条曲线围几次，而是所有特征轨迹围点的净圈数之和等于右半平面极点数。为什么是“和”？因为 (\det(I+L)=\prod_i(1+\lambda_i))，行列式的辐角是各因子辐角之和，辐角原理作用在乘积上，自然变成圈数相加。这个推导链条是理解广义奈奎斯特的关键：行列式判据和特征轨迹判据在数学上等价，但特征轨迹把“乘积抵消”的风险拆开了，每条轨迹单独可见，工程上更容易定位是哪一对特征值在逼近临界点。
这里有个容易忽略的细节：特征值在复平面上是连续变化的，但特征值的“编号”在频率扫描过程中可能发生交换。当两条特征轨迹在某个频率附近靠得很近时，数值计算出的特征向量会剧烈旋转，导致特征值排序跳变，画出来的曲线出现“断裂”或“交叉错位”。这不是系统本身不连续，而是特征分解的数值不稳定性。常见的处理办法是对特征向量做相位对齐，或者用更稳定的广义特征值求解器，保证轨迹在频率轴上连续跟踪。
2.3 为什么工程上更关心“离 -1 有多远”而不是“围没围”
判据只回答稳定与否，但实际调试时，工程师真正想知道的是“还有多少余量”。SISO 里用幅值裕度和相位裕度描述离 -1 点的距离，MIMO 里对应的概念是最小奇异值 (\sigma_{\min}(I+L(j\omega)))。当 (\sigma_{\min}) 接近 0 时，说明 (I+L) 接近奇异，闭环接近不稳定边界。广义奈奎斯特曲线和这个指标是互补的：曲线告诉你绕点关系，奇异值曲线告诉你最薄弱的方向和频率。很多翻车案例是判据显示稳定，但 (\sigma_{\min}) 在某个频段已经压到 0.1 以下，实际系统一加扰动就振荡。所以画广义奈奎斯特曲线时，我一般会同时把 (\sigma_{\min}(I+L)) 叠在同一张频率图上，两条信息一起看，才不会误判“稳定但脆弱”的状态。
3. 特征值轨迹怎么算出来：频率扫描、特征分解与连续性处理
3.1 频率扫描的离散化：步长、范围与对数分布
广义奈奎斯特曲线的数值实现，第一步是把连续频率轴离散化。设扫描范围为 ([\omega_{\min}, \omega_{\max}])，采样点数为 (N)。这里有两个参数直接决定曲线质量：步长和分布方式。线性等距采样在低频段分辨率过高、高频段分辨率不足，而奈奎斯特曲线在穿越 -1 点附近往往需要更细的采样才能看清绕向。常见的做法是对频率取对数分布：
% 频率扫描：对数分布，覆盖低频到高频 wmin = 1e-2; % 起始频率，rad/s wmax = 1e3; % 终止频率，rad/s N = 2000; % 采样点数 w = logspace(log10(wmin), log10(wmax), N); % 对数分布频率向量对数分布的好处是低频段点密、高频段点疏，符合大多数被控对象“低频增益变化快、高频趋于衰减”的特性。但要注意，如果系统在某个中频段有尖锐谐振，对数分布可能漏掉峰值附近的细节，这时需要在谐振频率附近局部加密。参数选择上，(\omega_{\min}) 一般取闭环带宽的 1/100 以下，(\omega_{\max}) 取开环增益降到 -40dB 以下的频率，保证曲线首尾都落在远离 -1 点的区域，绕向统计才完整。
3.2 逐频点特征分解：eig 的调用与复数矩阵处理
对每个频率点 (\omega_k)，计算 (L(j\omega_k)) 的特征值。MATLAB 里直接用eig：
% 假设 L 为 m×m 的 LTI 系统，先转成频率响应矩阵 m = size(L, 1); % 输入输出维数 Lf = squeeze(freqresp(L, w)); % m×m×N 的频率响应数组 Lf = reshape(Lf, m, m, []); % 统一成三维数组，兼容单频点情况 lambda = zeros(m, N); % 预分配特征值存储 for k = 1:N Lk = Lf(:, :, k); % 第 k 个频点的复矩阵 ek = eig(Lk); % 特征值分解，返回 m 个复特征值 lambda(:, k) = ek; % 按列存入 end这段代码的关键在freqresp的输出维度处理。当系统是单输入单输出时，freqresp返回的是 (1\times1\times N)，squeeze之后变成 (N\times1)，直接reshape成 (m\times m\times N) 会出错。所以要先squeeze再按 (m) 重塑，保证维度一致。另一个坑是eig对复数矩阵返回的特征值顺序不保证稳定，相邻频点之间可能交换位置，这就是 2.2 节提到的轨迹跳变问题。解决办法是在循环里做排序对齐：以上一个频点的特征值为参考，用最小距离匹配当前频点的特征值，保证第 (i) 条轨迹在频率轴上连续。
3.3 轨迹连续性：特征值排序与相位展开
排序对齐的常见实现是匈牙利算法或贪心最近邻匹配。贪心版本足够用：
% 特征值轨迹连续性对齐：贪心最近邻 lambda_sorted = zeros(m, N); lambda_sorted(:, 1) = sort(lambda(:, 1)); % 首点按实部排序 for k = 2:N prev = lambda_sorted(:, k-1); % 上一频点已对齐的特征值 cur = lambda(:, k); % 当前频点原始特征值 used = false(m, 1); % 标记已匹配 for i = 1:m d = abs(cur - prev(i)); % 到上一轨迹的距离 d(used) = inf; % 排除已用过的特征值 [~, idx] = min(d); % 最近邻匹配 lambda_sorted(i, k) = cur(idx); used(idx) = true; end end这段逻辑的意思是：每条轨迹在相邻频点之间不应该“跳远”，距离最近的特征值大概率是同一条轨迹的延续。贪心匹配在特征值靠得很近时仍可能出错，但配合足够密的频率采样，误差会被限制在局部。更严谨的做法是对特征向量做相位对齐后再匹配，因为特征值的连续性本质上由特征向量的连续性决定。实际画图时，如果看到某条轨迹在某个频点突然折返，先怀疑排序问题，再怀疑系统本身有模态交换。
提示：特征值排序只影响曲线的“连线方式”，不影响绕点圈数的统计结果。但如果要逐条轨迹分析裕度，排序错误会让裕度归属到错误的模态上，排查时容易误判。
4. 绕点圈数怎么数：从曲线到稳定性结论的最后一公里
4.1 包围圈数的数值统计：辐角增量法
有了连续的特征轨迹，下一步是统计每条轨迹逆时针包围 (-1+j0) 的净圈数。最直接的方法是计算轨迹相对 -1 点的辐角总增量：
[ N_i = \frac{1}{2\pi} \Delta \arg\left(\lambda_i(j\omega) - (-1)\right) ]
离散实现时，对相邻频点的辐角差做累加，并处理 (-\pi) 到 (\pi) 的跳变：
% 统计第 i 条轨迹包围 -1 点的净圈数 zi = lambda_sorted(i, :) - (-1); % 平移到以 -1 为原点 ang = angle(zi); % 各频点辐角，范围 (-pi, pi] dang = diff(ang); % 相邻辐角差 dang(dang > pi) = dang(dang > pi) - 2*pi; % 处理正跳变 dang(dang < -pi) = dang(dang < -pi) + 2*pi; % 处理负跳变 Ni = sum(dang) / (2*pi); % 净圈数，逆时针为正这里的关键是跳变修正。angle返回的辐角被折叠在 ((-\pi,\pi])，相邻两点真实相位差如果超过 (\pi)，直接相减会得到错误的大跳变。修正逻辑是把超过 (\pi) 的差值减去 (2\pi)、低于 (-\pi) 的加上 (2\pi)，还原真实增量。这个细节在频率采样较疏时尤其重要，采样越疏，单步相位变化越大，漏修正就会把圈数算错。血泪经验是：先用较密的采样跑一遍，确认圈数稳定，再逐步放宽采样看结果是否漂移，漂移说明步长已经不够。
4.2 判据落地：P 的确定与稳定结论的表述
圈数算出来后，还要确定 (P)，即开环传函矩阵在右半平面的极点数。对有理传函矩阵，(P) 等于所有元素通分后分母多项式的右半平面根数，计重数。数值上可以用pole求极点再数实部大于零的个数，但要注意数值误差导致的“虚轴附近极点”误判。常见做法是设一个容差 (\epsilon)，实部大于 (\epsilon) 才算右半平面极点，(\epsilon) 取 (10^{-6}) 量级。最终结论写成：
[ \sum_i N_i = P \quad \Rightarrow \quad \text{闭环稳定} ]
如果左边不等于右边，差值就是闭环在右半平面的极点数。这里有个容易混淆的点：广义奈奎斯特判据统计的是“净圈数之和”，不是“每条轨迹各自等于 P”。某条轨迹可能顺时针围了 2 圈、另一条逆时针围了 3 圈，净和是 1，只要等于 P 就稳定。逐条看会误以为有轨迹“围错方向”就不稳，实际上方向本身不构成判据，净和才是。
4.3 和奇异值曲线的交叉验证
单靠圈数统计，在特征轨迹贴近 -1 点时数值上很脆弱——轨迹稍微抖动，圈数就可能从 0 跳到 1。所以工程上我习惯用 (\sigma_{\min}(I+L(j\omega))) 做交叉验证：
% 最小奇异值曲线，作为稳定裕度的辅助判据 sv = zeros(N, 1); for k = 1:N Ik = eye(m) + Lf(:, :, k); % I + L(jw) sv(k) = min(svd(Ik)); % 最小奇异值 end semilogx(w, 20*log10(sv)); % dB 显示当 (\sigma_{\min}) 在某个频段明显下凹但未到 0，说明系统稳定但裕度薄；如果 (\sigma_{\min}) 穿越 0dB 以下，基本可以判定闭环不稳或临界稳定，此时再回头看特征轨迹，通常能找到某条轨迹正好穿过 -1 点附近。两条信息互相印证，比单看圈数可靠得多。参数上，(\sigma_{\min}) 的谷值频率往往就是系统最脆弱的频段，调控制器时优先在这个频段加阻尼或降增益。
5. 避坑与排查：广义奈奎斯特绘制里最容易翻车的几件事
5.1 现象：曲线画出来是散的，轨迹在某个频点突然断裂
原因：eig返回的特征值顺序在相邻频点间发生交换，排序对齐没做或做错了。特征向量在模态接近时剧烈旋转，导致特征值编号跳变。解决：按 3.3 节的贪心最近邻做轨迹对齐，并检查频率采样是否过疏；如果对齐后仍断裂，把该频段局部加密到原来的 5 到 10 倍再观察。
5.2 现象：圈数统计结果随采样点数变化，加密后从 0 变成 1
原因：轨迹在 -1 点附近有小的“毛刺”或数值抖动，稀疏采样时毛刺被跨过，密集采样时毛刺被捕捉，辐角增量法把毛刺算成了额外一圈。解决：先确认轨迹是否真的穿过 -1 点附近，用 (\sigma_{\min}) 曲线看该频点是否接近奇异；如果是数值毛刺，对特征值做轻度平滑或提高eig的数值精度（避免在病态矩阵上直接分解，可先做平衡化）。
5.3 现象：判据显示稳定，但实际闭环一加扰动就振荡
原因：只看了圈数，没看裕度。特征轨迹虽然没围住 -1，但某条轨迹离 -1 很近，(\sigma_{\min}(I+L)) 在该频段已经很低，系统处于“稳定但脆弱”的状态。解决：把 (\sigma_{\min}) 曲线和特征轨迹叠画，检查谷值频率和谷值大小；谷值低于 0.2 时就要考虑加控制器阻尼或降低该频段增益。
5.4 现象：P 算出来是 0，但系统明明有右半平面极点
原因：pole求出的极点实部因数值误差落在虚轴附近，容差设得太小或太大都会误判。解决：设 (\epsilon=10^{-6}) 量级的容差，并交叉检查分母多项式的根；对接近虚轴的极点，单独打印实部数值人工确认，不要直接信sum(real(p)>0)的布尔统计。
5.5 现象：MIMO 系统维数稍高，freqresp输出维度对不上，reshape报错
原因：单频点或单输入单输出时freqresp的返回维度退化，squeeze后维度数变化，直接reshape成 (m\times m\times N) 会错位。解决：先判断size再重塑，或用permute把频率维统一放到最后一维；对 (m=1) 的退化情况单独走标量分支，不要硬套矩阵流程。
6. 设计取舍：为什么特征轨迹比行列式轨迹更值得画
回到最初的问题：广义奈奎斯特曲线到底“广义”在哪。表面看是把一条曲线变成多条，深层看是把一个乘积判据拆成了可归因的因子判据。行列式轨迹 (\det(I+L(j\omega))) 在数学上完全等价，画起来也只有一条线，看起来更省事。但它的问题在于“乘积会掩盖个体”：两个特征值分别从 -1 的两侧靠近，乘积可能始终远离原点，行列式轨迹显示一切正常，而实际系统已经处在模态临界。特征轨迹把每个因子单独摊开，哪条轨迹在逼近 -1、逼近到什么程度，一目了然。这个取舍的代价是计算量和绘图复杂度上升，还要处理特征值排序的数值问题，但换来的是可归因性——出问题时你知道该调哪个模态，而不是对着一张“看起来没问题”的行列式曲线干瞪眼。
另一个取舍是判据和裕度的关系。圈数判据是二值的，稳定或不稳定，但工程上真正有价值的是“离不稳定还有多远”。所以我一般不会只画特征轨迹就收工，而是把 (\sigma_{\min}(I+L)) 作为固定搭配一起看。这两者的关系有点像体检报告里的“有没有病”和“各项指标离正常范围边界多远”，前者给结论，后者给行动依据。特征轨迹负责回答“围没围住”，奇异值曲线负责回答“最薄弱的频段和方向在哪”，缺一个都容易在调试时走弯路。
验证方法上，我习惯用一个小规模已知系统做回归：构造一个 2×2 的开环矩阵，故意让其中一个特征值在某个频率穿过 -1 附近，看绘制的轨迹和圈数统计是否符合预期，再和 (\sigma_{\min}) 的谷值位置对照。这个回归用例不需要复杂，但每次改完排序逻辑或频率采样策略都跑一遍，能挡住大部分数值层面的翻车。从那以后我每次画广义奈奎斯特曲线，都强制先把特征值连续性对齐和 (\sigma_{\min}) 交叉验证走一遍，再去看圈数结论——顺序反了，很容易被一条“看起来围住了”的曲线带偏。希望帮到你。
本文还有配套的精品资源，点击获取
