For those who don't have Control Toolbox, let's see how to draw a Bode plot with only basic Matlab functions.
Assuming that a transfer function is as below.
1
-------
s + 1
Then, the magnitude and phase of the transfer function are neccsary for bode plot.
The magnitude and phase can be calculated by replacing s with jw.
where, j is the imaginary unit of complex number, and w is the frequency.
1 1
------- --> -----------
s + 1 j*w + 1
--------------------------------------------------------
% Define frequency
freq = 0.01:0.01:1000;
% Transfer function
tf = 1 ./ (1 + i*freq);
% Magnitude
m = 20 * log10(abs(tf));
% Phase
phase = angle(tf);
% Plot
subplot(2,1,1)
semilogx(freq,m);
grid on
ylabel('Magnitude (dB)');
subplot(2,1,2)
semilogx(freq,phase);
grid on
ylabel('Phase (deg)');
xlabel('Frequency (rad/sec)');
--------------------------------------
As a result of the code above,
-----------------------------------------------------------------------
General method for Bode plot
% Define the numerator and denominator
num=[-0.1, -2.4, -181, -1950];
den = [1, 3.3, 990,2600];
tnum = [];
tden = [];
an = length(num);
bn = length(den);
% Numerator
for N = 1:1:an
if num(N) > 0
si ='+';
elseif num(N) <0 p="">
si = '';
end
temp = sprintf('%s%d%s%s%d ',si,num(N),'*','s.^',an-N);
tnum = [tnum temp];
end
% Denominator
for N = 1:1:length(den)
if den(N) > 0
si ='+';
elseif den(N) <0 p="">
si = '-';
end
temp = sprintf('%s%d%s%s%d ',si,den(N),'*','s.^',bn-N);
tden = [tden temp];
end
% Define a system transfer function
sys = @(s) eval(tnum)./ eval(tden);
% Define frequency as x
w=0.01:0.01:1000;
% Magnitude (dB)
amp=20*log10(abs(sys(w*1i)));
% Phase
ang=angle(sys(w*1i))*180/pi;
% Plot using subplot
subplot(211)
semilogx(w,amp);
ylabel('Magnitude (dB)')
grid on
subplot(212)
semilogx(w,ang);
xlabel('frequency (rad/sec)')
ylabel('phase (deg)')
grid on
