matlab画的bode图能看看频率从-20dB/dec变为-40dB/dec的交接频率吗？怎么看
- 提问者网友：爱了却不能说
- 2021-01-10 15:42
- 二级知识专家网友：低音帝王
- 2021-01-10 16:30
有些情况下，可以从两端较为平直的部分按照直线延伸得到交点，但如果模型中各转折频率相距比较近的话，幅相曲线会很难看出平直的部分，这种情况下上述方法不适用。
可以通过求零极点的方法直接确定各典型环节的转折频率，并按照大小排序，进而得到各转折频率及对应的斜率变化。这里提供一个我多年前编写的画折线图的函数，供参考。
效果图：
参考代码：
function corner( sys )
% Corner plot
if ~nargin, sys = zpk([-2 -5], [-1+i -1-i 0 -3], 1); end
KeyFun = [ ...
'AsciiVal= abs(get(gcf,''CurrentCharacter''));' ...
'if ~isempty(AsciiVal),' ...
' if AsciiVal==27,' ...
' close(gcf);' ...
' end,' ...
'end' ];
if 0,
bode(sys)
ax=findobj(gcf,'type','axes');lines=findobj(gcf,'type','line');set(lines,'linewidth',1.5)
set(ax,'units','points');set(gcf,'units','points');set(gcf,'pos',get(gcf,'pos')-[0 0 0 30])
set(gcf, 'KeyPressFcn', KeyFun )
hold on;
end
sys
[mag, phase, w] = bode(sys);
[p, z] = pzmap( sys );
[z, p, K] = zpkdata( sys );
z = z{1};
p = p{1};
K = K * prod( -z(z~=0) ) / prod( -p(p~=0) );
v = sum( p == 0 ) - sum( z == 0 );
p = p( p ~= 0 );
z = z( z ~= 0 );
wp = ones(0,2);
for i = 1 : length(p),
if imag(p(i)) == 0,
% for non-zero real poles
wp = [wp; abs(p(i)) -20];
elseif imag(p(i)) > 0,
% for complex conjugate poles, only consider those have positive imaginary parts
wp = [wp; abs(p(i)) -40];
end
end
wz = ones(0,2);
for i = 1 : length(z),
wz = [wz; abs(z(i)) 20];
end
wmin = w(1);
wmax = w(end);
wmin = min([wp(:,1); wz(:,1)]);
wmax = max([wp(:,1); wz(:,1)]);
if wmax == wmin, fac = 10; else, fac = 3; end
wmin = wmin / fac; wmax = wmax * fac;
ww = [wp; wz; wmin -20*v; wmax 0];
[mag, phase, w] = bode(sys, {wmin wmax});
[temp, idx] = sort( ww(:, 1) );
ww = ww(idx, :);
ww(:, 2) = cumsum( ww(:, 2) );
h(1) = 20 * log10(K) - 20 * v * log10(ww(1,1));
for i = 2 : size(ww, 1)
h(i) = h(i-1) + ww(i-1, 2) * log10( ww(i,1) / ww(i-1,1) );
end
fprintf(' v = %i, K = %.3g, 20lgK = %.3g
', v, K, 20 * log10(K) );
fprintf(' %i', ww(1, 2) );
for i = 2 : size(ww, 1) - 1,
fprintf(' --(%.3g)--> %i', ww(i, 1), ww(i, 2) );
end
fprintf('
');fig_h = 400;
fig_v = 300;
ax_l = 60;
ax_r = 20;
ax_b = 40;
ax_off = 15;
ax_t = 15;
ax_h = fig_h - ax_l - ax_r;
ax_v = ( fig_v - ax_b - ax_t - ax_off ) / 2;
figure('Name', 'Corner plot', 'KeyPressFcn', KeyFun, 'Units', 'Points', 'Pos', [100 100 fig_h fig_v], 'Color', 'w', 'Number', 'off' );
%subplot(2, 1, 1);
axes( 'Units', 'Points', 'Pos', [ax_l fig_v-ax_v-ax_t ax_h ax_v] );
plot([ww(1,1) ww(end,1)], [0 0], 'Color', ones(1,3)*.8);
hold on
plot(ww(:,1), h(:), 'LineWidth', 2, 'color', 'k');
plot(w, 20*log10(permute(mag, [3 1 2])), ':', 'color', 'k');
for i = 2 : length(h) - 1,
plot([ww(i,1) ww(i,1)], [h(i) 0], ':');
if h(i) >= 0,
valign = 'top';
halign = 'right';
else
valign = 'bottom';
halign = 'left';
end
text(ww(i,1), 0, num2str(ww(i,1), '%.3g'), 'VerticalAlignment', valign, 'HorizontalAlignment', halign);
end
set(gca, 'XScale', 'log', 'XtickLabel', '', 'Xlim', [wmin wmax], 'Box', 'on');
% convert 'Points' to 'Data'
p2d_v = ax_v / diff(get(gca, 'ylim'));
p2d_h = ax_h / diff(log10(get(gca, 'xlim')));
for i = 1 : length(h) - 1,
x = sqrt(ww(i,1)*ww(i+1,1));
y = (h(i)+h(i+1))/2;
ang = atan2( (h(i+1) - h(i)) * p2d_v, log10( ww(i+1,1)/ww(i,1) ) * p2d_h ) * 180 / pi;
if y >= 0,
valign = 'bottom';
y = y + 2 / p2d_v;
else
valign = 'top';
y = y - 1 / p2d_v;
end
text(x, y, sprintf('[%.3g]', ww(i,2)), 'VerticalAlignment', valign, 'HorizontalAlignment', 'center', 'Rotation', ang);
end
ylabel('Magnitude (dB)');
%subplot(2, 1, 2);
axes( 'Units', 'Points', 'Pos', [ax_l ax_b ax_h ax_v] );
%plot([ww(1,1) ww(end,1)], [-180 -180], 'Color', ones(1,3)*.8);
hold on
plot(w, permute(phase, [3 1 2]), 'LineWidth', 2, 'color', 'k');
set(gca, 'XScale', 'log', 'Xlim', [wmin wmax], 'Box', 'on');
ylabel('Phase (deg)');
xlabel('Frequency (rad/sec)');
set( findobj(gcf, 'type', 'axes'), 'Units', 'normal' )
print -dmeta
追问赞追答十多年前编的，现在看来有不少不足之处，但基本的意思应该是有了。
- 1楼网友：执傲
- 2021-01-10 18:08