Monday, August 14, 2017

[Matlab] Bode plot without Control Toolbox

When it comes to Bode plot, it is easy to draw a Bode plot with control toolbox, but Not everybody can get this toolbox.
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

Possible matlab code for bode plot is as below.
--------------------------------------------------------
% 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