There are are a few output variabales missing in the fcn: bode_as21. Besides, Wl is incorrectly assigned that should be tf and called by the variable name of SYS. Here is the corrected code:
clc
close all
clear all
% Початкові дані:
k=10;
T1=1;
T2=0.5;
% Вектор абцис частотних функцій:
omega=linspace(0,5);
p=omega*1i;
% Передатна функція:
W=k./(T2*T1*p.^2+T1*p+1);
% ДЧХ:
Wre=real(W);
% УЧХ:
Wim=imag(W);
% Модуль ПФ:
A=abs(W);
% Кут (і перевод в градуси):
fi=angle(W);
fi_d=fi.*180./pi;
% a)ДЧХ і УЧХ в одному вікні;
plot(omega, Wre, omega, Wim), grid on;
xlabel('ω');
ylabel('Дійсна та уявна частини');
title('ДЧХ та УЧХ');
legend('ДЧХ','УЧХ');
% б) АЧХ і ФЧХ у двох підвікнах ...
...одного вікна (одне під іншим);
figure(2);
subplot(2,1,1);
plot(omega, A), grid on;%АЧХ;
title('АЧХ');
xlabel('ω');
ylabel('A');
subplot(2,1,2);
plot(omega, fi_d), grid on;%ФЧХ;
title('ФЧХ');
xlabel('ω');
ylabel('φ');
% в) Bode-діаграми (ЛАЧХ і ЛФЧХ) у двох підвікнах одного вікна (одне...
...під іншим):
% Вектор абцис частотних функцій і кут:
lomega=logspace(-1,1);
pl=lomega*1i;
% Початкові дані:
k=10;
T1=1;
T2=0.5;
p=pl;
% Передатна функція:
Wl=k./(T2*T1*p.^2+T1*p+1);
% Модуль передатної функції і її логарифмічний вигляд:
Al=abs(Wl);
L=20*log10(Al);
% Перевод даних в градуси:
fil=angle(Wl);
fil_d=fil.*180./pi;
H=20*log10(fil_d);
figure(3);
subplot(2,1,1);
semilogx (lomega, L), grid on
title('ЛАЧХ');%ЛАЧХ;
xlabel('ω');
ylabel('A');
subplot(2,1,2);
semilogx (lomega, H), grid on
title('ЛФЧХ');%ЛФЧХ;
xlabel('ω');
ylabel('φ');
figure(4);
obj=tf([k],[T2*T1 T1 1]);
bode(obj), grid on;%Bode-діаграми методом функції bode;
%г) діаграму Найквіста (АФЧХ) у полярних координатах, у Декартових...
...координатах та за допомогою функції nyquist:
figure(5);
polar(fi, A), grid on;%у полярних координатах;
title('АФЧХ');
xlabel('A');
ylabel('φ');
figure(6);
plot(Wre,Wim), grid on;%у Декартових координатах;
SYS=tf(k,[T2*T1 T1 1]);
title('АФЧХ');
xlabel('Дійсна частина');
ylabel('Уявна частина');
figure(7);
nyquist(SYS);%За допомогою функції nyquist;
% На ЛАЧХ нанести асимптотичну характеристику (пізніше):
[wt,wlt,knt,Lt,wc] = bode_as21(SYS, 1e-4, 1e4)
function [wt, wcl, knt, La, wc] = bode_as21 (sys, w0, wk)
[Z, P, K] = zpkdata(sys, 'v'); % Визначення нулів та полюсів
wp = abs (P); % Амплітуди (частоти) полюсів
kp = length (P); n = kp; % Кількість полюсів
wz = abs(Z); % Амплітуди (частоти) нулів
kz = length(Z); m = kz; % Кількість нулів
% Вилучення із масиву частот нульових (нейтральних) полюсів та нулів:
ip0 = find (wp == 0); % Порядкові номери нейтральних нулів та полюсів
iz0 = find (wz == 0); % Порядкові номери нейтральних нулів та полюсів
wp (ip0) = []; wz (iz0) = [];
kp0 = length (ip0); kz0 = length (iz0); n1 = n - kp0; m1 = m - kz0;
kn0 = (kz0 - kp0) * 20; % Кількість нейтральних нулів та полюсів
w = [wz', wp']; % Об’єднання частот спряження в один вектор-рядок
dkn = [ones(1, kz-kz0), -ones(1, kp-kp0)]*20; % Визначення зміни нахилів ЛАЧХ
% кожній частоті
[ws, is] = sort (w); % Сортування частот спряження за зростанням
dkns = dkn (is); % Перестановка змін нахилів відповідно до перестановки
% частот
kns = kn0 + cumsum (dkns); % Накопичення нахилів ЛАЧХ зі врахуванням
% нахилу початкової ділянки
% Об’єднання кратних частот:
ik = find (diff (ws) ==0);
ws (ik) = [];
kns (ik) = [];
kf =K* prod (wz) / prod (wp); % Фіктивний коефіцієнт підсилення
wf = kf ^ (-20/kn0); % Фіктивна частота зрізу
% Визначення діапазона частот за замовчанням
if nargin ==1 % Якщо діапазон частот не задано користувачем при зверненні
% до функції
wsl = log10(ws); % Логарифми частот спряження
if kn0==0, wl_min = wsl(1); else wfl = log10(wf); wl_min = min (wfl, wsl(1)); end
wl0 = floor(wl_min); wlk=ceil(wsl(end)); % Логарифми меж діапазону частот
w0=10^wl0; wk=10^wlk; % Діапазон частот
end
if kn0 == 0
L0 = 20 * log10 (kf);
else
L0 = - kn0 * log10 (wf / w0);
end
% Формування таблиці для побудови АЛАЧХ:
wt = [w0 ws wk]; % Приєднання до частот спряження граничних
% частот
wtl = log10 (wt); % Логарифми об’єднаного вектору частот
d_wl = diff (wtl); % Різниці між логарифмами сусідніх частот
knt = [kn0 kns kns(end)]; % Приєднання до вектору коефіцієнтів нахилу
% початкового та кінцевого елементів
% Розрахунок ординат АЛАЧХ:
kt = length (wt); La = zeros (1, kt); La(1) = L0;
for i = 2 : kt
La(i) = La(i-1) + knt (i-1) * d_wl (i-1);
end
ic = find (La (1 : kt-1) .* La (2 : kt) < 0);
% Пошук частоти зрізу методом зворотного інтерполювання
wcl = interp1 (fliplr (La(ic : ic+1)), fliplr (wtl (ic : ic+1)), 0);
wc=10^wcl;
[A, fi, wl] = bode (sys,{w0,wk}); % Розрахунок амплітуд, фаз та частот для
% побудови точних частотних характеристик
A = A(:); fi = fi(:); % Позбавлення від зайвих індексів
L = 20*log10(A); % Розрахунок амплітуди в дБ
figure, semilogx (wl, L, wt, La, '--', 'LineWidth', 1.5 ), grid on, hold on % Точна ЛАЧХ
% та АЛАЧХ
semilogx ([w0,wk], [0 0], 'k') % Вісь частот
for i = 2 : kt-1
semilogx ([wt(i) wt(i)], [0 La(i)], 'r--', 'LineWidth',1.5) % Проекції точок спряження
% на вісь частот
end
Lmax = max (L); Lmin = min(L);
ylim ([Lmin-10 Lmax+10])
xlabel ('lg (omega)')
ylabel ('20lg |W (j*omega)|')
legend ('ЛАЧХ', 'АЛАЧХ')
end
