MATLAB 系统仿真与建模（十）—— 频域响应分析

日期：2020.12.30

参考：

1. 控制理论 MATLAB 教程

1 绘制伯德图（bode）

bode 图是直角坐标图

bode 命令可以计算连续线性时不变（定常）系统频域响应的幅值和相角，示例代码如下：

bode(num,den) 
bode(num,den,w) 
bode(A,B,C,D)
bode(A,B,C,D,w)
bode(A,B,C,D,iu,w) 
bode(sys)

不带左端参数的时候屏幕上会直接绘制出伯德图。图中幅值单位为分贝（dB）

当使用左端参数：

[mag, phase, w] = bode(num, den, w); % 其他同理

矩阵mag 和 phase 分别包含系统频率响应在指定频率点上算得的幅值和相角。相角的单位为度 (°) 。 幅值可以通过如下语句转换为分贝值：

magdB = 20*log10(mag);

如果需要指定频率范围，可使用命令：

w = logspace(d1, d2); % 在十倍频程 10^d1 和 10^d2 内生成一个包含50个点的，呈对数均匀分布的向量（包含起点和终点）
w = logspace(d1, d2, n); % 十倍频程 10^d1 和 10^d2 内 生成呈对数均匀分布的n个点(n个频率点包含起点和终点）

eg. 如下传递函数，绘制 bode 图：
$$G(s)=\frac{25}{s^2+4s+25}$$
代码如下：

num = [0 0 25];
den = [1 4 25];
G = tf(num, den);
bode(G)
grid
title('Bode Diagram of G(s)=25/(s^2+4s+25)')

绘制出图片如下：

如果期望幅值的范围介于某一范围之间（例如-50dB~50dB之间），就可以在伯德图上添加幅值为-50dB和50dB两条直线。

如果期望相角的范围在-150°~150°之间，举个例子：

系统如下：
$$G(s)=\frac{9(s^2+0.2s+1)}{s(s^2+1.2s+9)}$$
我们规定范围内绘制bode图，代码如下：

num = [9 1.8 9];
den = [1 1.2 9 0];
w = logspace(-2, 3, 100);
[mag, phase, w] = bode(num, den, w);
magdB = 20*log10(mag);
dBmax = 50*ones(1,100); 
dBmin = -50*ones(1,100); %点数必须与 w 中频率点的数目相等

% ********************************************************************
subplot(2, 1, 1)
semilogx(w, magdB, 'o', w, magdB, '-', w, dBmax, '--', w, dBmin, ':')、
grid
title('BodeDiagram of G(s)=9(s^2+0.2s+1)/(s(s^2+1.2s+9)]')
xlabel('Frequency (rad/sec)')
ylabel('Gain dB')

% ********************************************************************
subplot(2, 1, 2)
pmax = 150*ones(1,100); 
pmin = -150*ones(1,100);
semilogx( w, phase, 'o', w, phase, '-', w, pmax, '--', w, pmin, ':')
grid
xlabel('Frequency (rad/sec)')
ylabel('Phase deg')

2 绘制奈奎斯特图（nyquist）

nyquist 图是极坐标图

MATLAB 命令 nyquist 可以计算线性时不变连续系统的频域响应。如果调用时不带左端参数， nyquist 就在屏幕上绘制奈奎斯特图

代码如下：

nyquist(num, den);
nyquist(num,den,w); % 包含了用户指定的频率向量 w 
nyquist(A,B,C,D);
nyquist(A,B,C,D,w);
nyquist(A,B,C,D,iu,w);
nyquist(sys);

使用了上述任一命令，接下来再使用 grid 命令，则不会产生 x-y 网格线（水平和垂直方向的网格线）。如果需要的话可以使用带左端参数的 nyquist 指令，再配合使用 grid 命令来实现。

[re,im,w) = nyquist(num,den);
[re,im,w) = nyquist(num,den,w); 
[re,im,w) = nyquist(A,B,C,D);
[re,im,w) = nyquist(A,B,C,D,w);
[re,im,w) = nyquist(A,B,C,D,iu,w); 
[re,im,w) = nyquist(sys);

MATLAB 将系统的频域响应的数值赋予矩阵 re、im 和 w，但不在屏幕上绘制奈奎斯特图。矩阵 re 和 im 分别为系统在由向量 w 给定的频率点上的频域响应的实部和虚部。应当注意的是，re 和 im 的列数与系统输出的数目相同， 其每一行对应于 w 的一个元素。若要绘制奈奎斯特图， 则可以使用 plot 指令。

举个例子：
$$G(s)=\frac{1}{s^2+0.8s+1}$$
代码如下：

num = [1];
den = [1 0.8 1];
nyquist(num, den)
title('Nyquist Plot of G(s) = 1/(s^2+0.8s+1)')

绘制出图片如下：

如果希望依据自行设定的实轴和虚轴范围绘制奈奎斯特图，例如实轴和虚轴的范围均为从 -2 到 +2，就可以输入命令：

v = [-2 2 -2 2];
axis(v);

注意： 在绘制奈奎斯特图时，如果 MATLAB 运算涉及 “Divided by zero” （被零除），所绘制的奈奎斯特图就可能是错误的，或者不适合进一步分析的需要。

例如：
$$G(s)=\frac{1}{s(s+1)}$$
系统正常绘制出来的图为：

很明显这不利于分析，于是我们限制了坐标范围后 axis([-2 2 -5 5]) 绘制出的图片如下：

添加上 grid 则可以叠印等幅值圆：

3 绘制对数幅相特性图 / 尼柯尔斯图（nichols）

• 尼柯尔斯图（英语：Nichols plot）是将线性非时变系统在不同频率下的增益分贝值及相位绘在一直角坐标系的图上，尼柯尔斯图将二种波德图（波德增益图及波德相位图）结合成一张图，而频率只是曲线中的参数，不直接在图中显示。以系统的对数分贝幅值为纵坐标，以相角或相角裕度为横坐标绘制而成的图形。尼柯尔斯图的命名是来自美国控制工程师尼柯尔斯（Nathaniel B. Nichols）。

• 在伯德图中，$G(jw)$ 的频率响应特性在半对数坐标纸上用两条独立的曲线来表示，即对数幅值曲线和相角曲线。但是，在对数幅相特性图中，伯德图中的这两条曲线合并为一条曲线。通过从伯德图中读取对数幅值和相角的数值，可以很轻松地手工绘制对数幅相特性图。值得注意的是，在对数幅相特性图中绘制 $G(jw)$ 的增益常数变化时仅仅是将曲线向上（增益增加时）或向下（增益减小时）平移，但保持曲线的形状不变。

• 对数幅相特性图的优势在于，它可以快速判断闭环系统的相对稳定性，并很容易设计补偿器。

• 正弦传递函数 $G(jw)$ 和 1 / G ( j w ) 1/ G(jw)