Frequency Response & Bode Diagrams
Introduction
Consider an asymptotically stable LTI continuous-time that is fully described (that is, no common factors, i.e. no hidden dynamics) by the transfer function 
Moreover, consider 
The regime sinusoidal response (approx. equal to
for
, where
is the settling time of the step response of the system
) is: The frequency response of the system is defined as: 
With a Little Help from the 'help' Command!
Here is the documentation related to the Control System Toolbox function we will use in this live script:
help evalfr
--- help for DynamicSystem/evalfr ---
DynamicSystem/evalfr - Evaluate system response at specific frequency
evalfr is a simplified version of freqresp meant for quick evaluation of
the system response at any point in the complex plane.
Syntax
frsp = evalfr(sys,x)
Input Arguments
sys - Dynamic system
dynamic system model | model array
x - Point in complex plane
complex scalar
Output Arguments
frsp - Frequency response
complex scalar | complex array
Examples
Evaluate Discrete-Time Transfer Function
Evaluate Frequency Response of Identified Model at Given Frequency
Frequency Response of MIMO State-Space Model
See also bode, freqresp, sigma
Introduced in Control System Toolbox before R2006a
Documentation for DynamicSystem/evalfr
help freqresp
freqresp - Evaluate system response over a grid of frequencies
Use freqresp to evaluate the system response over a grid of frequencies.
Syntax
[H,wout] = freqresp(sys)
H = freqresp(sys,w)
H = freqresp(sys,w,units)
[H,wout,covH] = freqresp(sys,___)
Input Arguments
sys - Dynamic system
dynamic system model | model array
w - Frequencies
vector
units - Units of the values in the input frequency vector
rad/TimeUnit (default) | 'cycles/TimeUnit' | 'rad/s' | ...
Output Arguments
H - Frequency response values
array
wout - Output frequencies corresponding to the frequency response
vector
covH - Covariance of the frequency response
array
Examples
Frequency Response of SISO System
Compute Frequency Response of System
Compute Frequency Response on Specified Frequency Grid
Compute Frequency Response and Associated Covariance
See also evalfr, bode, nyquist, nichols, sigma, interp,
spectrum (System Identification Toolbox)
Introduced in Control System Toolbox before R2006a
Documentation for freqresp
Other uses of freqresp
help bode
bode - Bode frequency response of dynamic system
This MATLAB function computes the frequency response of dynamic system
model sys and returns the magnitude and phase of the response at each
frequency in the vector wout.
Syntax
[mag,phase,wout] = bode(sys)
[mag,phase,wout] = bode(sys,w)
[mag,phase,wout,sdmag,sdphase] = bode(sys,w)
bode(___)
Input Arguments
sys - Dynamic system
dynamic system model | model array
w - Frequencies
{wmin,wmax} | vector | []
Output Arguments
mag - Magnitude of system response
3-D array
phase - Phase of system response
3-D array
wout - Frequencies
vector
sdmag - Standard deviation of magnitude
3-D array | []
sdphase - Standard deviation of phase
3-D array | []
Examples
Bode Plot of Dynamic System
Bode Plot at Specified Frequencies
Compare Bode Plots of Several Dynamic Systems
Bode Plot with Specified Line Attributes
Obtain Magnitude and Phase Data
Magnitude and Phase of MIMO System
Bode Plot of Identified Model
Obtain Magnitude and Phase Standard Deviation Data of Identified Model
Bode Plot of Model with Complex Coefficients
See also bodeplot, freqresp, nichols, nyquist, step
Introduced in Control System Toolbox before R2006a
Documentation for bode
Other uses of bode
Moreover, the following custom functions are available:
addpath(genpath("BodeDiagram"))
help asympt_bode
*** asympt_bode(SISOmodel, om_vector) ***
[MagData,PhaseData,out_om_v] = asympt_bode(SISOmodel, om_vector)
Asymptotic Bode diagram of a transfer function.%
The SISOmodel can be any of the available LTI system models,
i.e. a tf object (preferred), an ss object or a zpk object.
In the event that SISOmodel is not a tf object, it is converted to tf.
om_vector --> angular frequency values used to draw the graphs
(OPTIONAL - if not present, a suitable frequency range
is estimated on the basis of zeros and poles. An
appropriate number of values for the angular pulsation
in the newly determined range is then chosen).
If MagData and PhaseData vectors are not requested in the output,
asymptotic modulus and phase graphs are plotted in graphs
on a semi-logarithmic scale.
The array out_om_v contains the angular frequency values
used to compute the asymptotic diagram of the frequency response.
Last Update: 2024/05/24
----------------------------------------------
help drawBodediagrams
drawBodediagrams( )
Plotting Bode diagrams of the frequency response
(asymptotic and actual superimposed diagrams) of
continuous-time transfer functions.
INPUT variables:
Gs: transfer function, described as TF object; MANDATORY variable
omega_values: angular frequency values used to draw the graphs
(OPTIONAL - if not present, a suitable frequency range
is estimated on the basis of zeros and poles. An
appropriate number of values for the angular pulsation
in the newly determined range is then chosen).
asBode_color: colour of asymptotic graphs
(array: three values between 0 and 1) - OPTIONAL
Default value: [1, 0, 0] <--> red
asBode_line_width: line thickness for asymptotic graphs - OPTIONAL
Default value: 2
asBode_line_style: line style for asymptotic diagrams
(see 'LineStyle' in plot( ) ) - OPTIONAL
Default value: '-'
bode_color: colour of actual graphs
(array: three values between 0 and 1) - OPTIONAL
Default value: [0, 0, 1] <--> blue
bode_line_width: line thickness for actual graphs - OPTIONAL
Default value: 1.5
bode_line_style: line style for actual diagrams
(see 'LineStyle' in plot( ) ) - OPTIONAL
Default value: '-'
fig_handler: The figure handler of a figure already in use,
where Bode diagrams are drawn.
If it does not exist, a new figure is created.
comp_asympt_vs_actual_ph_diagrsFLAG:
a logic flag indicating whether the 360-degree jump
between the asymptotic phase diagram and the actual
phase diagram should be compensated (if present).
Default value: true
----------
OUTPUT variables (optional): the handlers of the couple of axes
corresponding to the two subplots.
A simple Example: Evaluation of the Frequency Response at a single angular Frequency
Then
What about using MATLAB?
s = tf('s'); % the transfer function building helper
Gs = 1/(1+s); % the transfer function of G(s)
sVAL = 1i*2; % j2 <--> the **complex frequency** corresponding to the input u(t)
% refer to the documentation of evalfr( ) for details
FreqResp_Gs_2 = evalfr(Gs, sVAL)
FreqResp_Gs_2 = 0.2000 - 0.4000i
% the frequency response G(j2): it's a complex value
% the magnitude of the frequency response at the angular frequency 2 rad/s
magG2 = abs(FreqResp_Gs_2)
% the phase angle of FreqResp_Gs_2 <--> arg G(j2) <-- NB rad
argG2_rad = angle(FreqResp_Gs_2)
% converting from rad to deg
argG2_deg = rad2deg(argG2_rad)
Bode Diagrams: asymptotic and effective Diagrams using MATLAB
Example 1: A Couple of Complex Poles, varying the Damping Factor
Consider the LTI system described by the transfer function:

xiVAL = [0.05, 0.2, 0.5, 0.7, 0.9]; % the damping factor values
Nxi = numel(xiVAL); % how many damping factors?
% let's preallocate a data-structure for storing the different transfer functions
% generate the transfer functions, with different damping factors
numGs = 100; % the coeffs of the polynomial at the numerator in G(s)
denGs = [ (1) (20*xiVAL(n)) 100];
% the coeffs of the polynomial at the denominator in G(s),
Gs_xi{n} = tf(numGs, denGs); % let's define the transfer function G(s)
% and store it in the data-structure
omVALS = logspace(-1,2, 1e4);
% We want to evaluate the frequency response at 10000
% angular pulsation values between 10^-1 and 10^2 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolors = [0 0.4470 0.7410; ...
0.9290 0.6940 0.1250; ...
0.4940 0.1840 0.5560; ...
0.4660 0.6740 0.1880; ...
0.3010 0.7450 0.9330; ...
0.8500 0.3250 0.0980; ...
0.6350 0.0780 0.1840]; % some different colors
% for the actual diagrams
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
[hax1, hax2] = drawBodediagrams(Gs_xi{n}, omVALS, asBcolor, 2.5, '-', ...
Bcolors(n, :), 1.5, '-', hf);
legend(hax1, 'asympt.', ['$\xi = ',num2str(xiVAL(1))], ...
'', ['\xi = ',num2str(xiVAL(2))], ...
'', ['\xi = ',num2str(xiVAL(3))], ...
'', ['\xi = ',num2str(xiVAL(4))], ...
'', ['\xi = ',num2str(xiVAL(5))], ...
legend(hax2, 'asympt.', ['\xi = ',num2str(xiVAL(1))], ...
'', ['\xi = ',num2str(xiVAL(2))], ...
'', ['\xi = ',num2str(xiVAL(3))], ...
'', ['\xi = ',num2str(xiVAL(4))], ...
'', ['\xi = ',num2str(xiVAL(5))], ...
Note: the asymptotic diagrams of the transfer functions, obtained by varying the damping coefficient, coincide as they should.
Example 2: a Minimum-Phase LTI System with a Pole at the Origin
Consider the LTI system described by the transfer function:
Analyzing the transfer function, you get

addpath(genpath('BodeDiagram'));
s = tf('s'); % the transfer function building helper
Gs = 100*(1+10*s)/(s * (1+2*s) * (1+0.4*s+s^2) )
Gs =
1000 s + 100
-----------------------------
2 s^4 + 1.8 s^3 + 2.4 s^2 + s
Continuous-time transfer function.
Model Properties
omVALS = logspace(-3,2, 1e4);
% We want to evaluate the frequency response at 10000
% angular pulsation values between 10^-2 and 10^1 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolor = [0 0.4470 0.7410]; % a different color for the actual diagram
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
[hax1, hax2] = drawBodediagrams(Gs, omVALS, asBcolor, 2.5, '-', ...
legend(hax1, 'asympt.', 'actual', ...
legend(hax2, 'asympt.', 'actual', ...
Example 3: a Minimum-Phase LTI System with a Pole at the Origin
Consider the LTI system described by the transfer function:
Analyzing the transfer function, you get

s = tf('s'); % the transfer function building helper
Gs = 100*(1+10*s)/(s * (1+s)^2 )
Gs =
1000 s + 100
---------------
s^3 + 2 s^2 + s
Continuous-time transfer function.
Model Properties
omVALS = logspace(-2,2, 5e4);
% We want to evaluate the frequency response at 50000
% angular pulsation values between 10^-2 and 10^2 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
% asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
% Bcolor = [0 0.4470 0.7410]; % a different color for the actual diagram
% hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
drawBodediagrams(Gs, omVALS)
% [hax1, hax2] = drawBodediagrams(Gs
% );, omVALS, asBcolor, 2.5, '-', ...
% legend(hax1, 'asympt.', 'actual', ...
% legend(hax2, 'asympt.', 'actual', ...
Example 4: an Unstable RHP Pole
Consider the LTI system described by the transfer function:
Analyzing the transfer function, you get

s = tf('s'); % the transfer function building helper
Gs = 0.1*s*(1+s)/((1+0.2*s) * (1 + 5*s)^2 *(1 - 0.1*s) )
Gs =
-0.1 s^2 - 0.1 s
------------------------------------------
0.5 s^4 - 2.3 s^3 - 25.98 s^2 - 10.1 s - 1
Continuous-time transfer function.
Model Properties
omVALS = logspace(-3,3, 5e4);
% We want to evaluate the frequency response at 50000
% angular pulsation values between 10^-2 and 10^2 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolor = [0 0.4470 0.7410]; % a different color for the actual diagram
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
[hax1, hax2] = drawBodediagrams(Gs, omVALS, asBcolor, 2.5, '-', ...
legend(hax1, 'asympt.', 'actual', ...
legend(hax2, 'asympt.', 'actual', ...
Example 5: Complex Poles with Different Damping Factor, but with the same Natural Frequency.
Consider the LTI systems described by the transfer functions:
Analyzing the transfer functions, you get

s = tf('s'); % the transfer function building helper
Gs1 = 50 * (1+0.4*s)/( (1+10*s) * (1 + 0.2*s +s^2) )
Gs1 =
20 s + 50
---------------------------
10 s^3 + 3 s^2 + 10.2 s + 1
Continuous-time transfer function.
Model Properties
Gs2 = 50 * (1+0.4*s)/( (1+10*s) * (1 + 1.6*s +s^2) )
Gs2 =
20 s + 50
----------------------------
10 s^3 + 17 s^2 + 11.6 s + 1
Continuous-time transfer function.
Model Properties
omVALS = logspace(-2,1, 1e4);
% We want to evaluate the frequency response at 10000
% angular pulsation values between 10^-2 and 10^1 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0]; % red color (RGB) for the asymptotic diagram
Bcolor = [0 0.4470 0.7410; ...
0.9290 0.6940 0.1250]; % a different color for the actual diagram
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
drawBodediagrams(Gs1, omVALS, asBcolor, 2.5, '-', ...
Bcolor(1,:), 1.5, '-', hf);
[hax1, hax2] = drawBodediagrams(Gs2, omVALS, asBcolor, 2.5, '-', ...
Bcolor(2,:), 1.5, '-', hf);
legend(hax1, 'asympt.', 'actual G^{''}(s)', ...
'', 'actual G^{"}(s)',...
legend(hax2, 'asympt.', 'actual G^{''}(s)',...
'', 'actual G^{"}(s)', ...
Note: the asymptotic diagrams of the frequency response of both the transfer functions coincide, as expected.
Example 6 - Non Minimum Phase System: a RHP Zero
Consider the LTI systems described by the transfer functions:
Analyzing the transfer functions, you get

s = tf('s'); % the transfer function building helper
Gs1 = (1+s)/( s^2+5*s+6)
Gs1 =
s + 1
-------------
s^2 + 5 s + 6
Continuous-time transfer function.
Model Properties
Gs2 = (1-s)/( s^2+5*s+6)
Gs2 =
-s + 1
-------------
s^2 + 5 s + 6
Continuous-time transfer function.
Model Properties
omVALS = logspace(-2,1, 1e4);
% We want to evaluate the frequency response at 10000
% angular pulsation values between 10^-2 and 10^1 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0; 0 1 0]; % red or green color (RGB) for the asymptotic diagrams
Bcolors = [0 0.4470 0.7410; ...
0.9290 0.6940 0.1250; ...
0.4940 0.1840 0.5560; ...
0.4660 0.6740 0.1880; ...
0.3010 0.7450 0.9330; ...
0.8500 0.3250 0.0980; ...
0.6350 0.0780 0.1840]; % a different color for the actual diagram
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
drawBodediagrams(Gs1, omVALS, asBcolor(1,:), 2.5, '-', ...
Bcolors(1,:), 1.5, '-', hf);
drawBodediagrams(Gs2, omVALS, asBcolor(2,:), 2.5, '-', ...
Bcolors(2,:), 2.5, '-', hf, true);
[hax1, hax2] = drawBodediagrams(Gs2, omVALS, asBcolor(2,:), 2.5, '-', ...
Bcolors(3,:), 1.5, '-', hf, false);
legend(hax1, 'asympt. G^{''}(s)', 'actual G^{''}(s)', ...
'asympt. G^{"}(s)', 'actual comp. G^{"}(s)',...
'','actual no-comp. G^{"}(s)','Location', 'best');
legend(hax2, 'asympt. G^{''}(s)', 'actual G^{''}(s)', ...
'asympt. G^{"}(s)', 'actual comp. G^{"}(s)',...
'','actual no-comp. G^{"}(s)','Location', 'best');
Note: the magnitude diagrams of the frequency response coincide, as expected.
Example 7 - Non Minimum Phase System: Negative Static Gain
Consider the LTI systems described by the transfer functions:
Analyzing the transfer functions, you get

s = tf('s'); % the transfer function building helper
Gs1 = (1+s)/( s^2+5*s+6)
Gs1 =
s + 1
-------------
s^2 + 5 s + 6
Continuous-time transfer function.
Model Properties
Gs2 = -(1+s)/( s^2+5*s+6)
Gs2 =
-s - 1
-------------
s^2 + 5 s + 6
Continuous-time transfer function.
Model Properties
omVALS = logspace(-2,1, 1e4);
% We want to evaluate the frequency response at 10000
% angular pulsation values between 10^-2 and 10^1 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0; 0 1 0]; % red or green color (RGB) for the asymptotic diagrams
Bcolors = [0 0.4470 0.7410; ...
0.9290 0.6940 0.1250; ...
0.4940 0.1840 0.5560; ...
0.4660 0.6740 0.1880; ...
0.3010 0.7450 0.9330; ...
0.8500 0.3250 0.0980; ...
0.6350 0.0780 0.1840]; % a different color for the actual diagram
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
drawBodediagrams(Gs1, omVALS, asBcolor(1,:), 2.5, '-', ...
Bcolors(1,:), 1.5, '-', hf);
drawBodediagrams(Gs2, omVALS, asBcolor(2,:), 2.5, '-', ...
Bcolors(2,:), 2.5, '-', hf, true);
[hax1, hax2] = drawBodediagrams(Gs2, omVALS, asBcolor(2,:), 2.5, '-', ...
Bcolors(3,:), 1.5, '-', hf, false);
legend(hax1, 'asympt. G^{''}(s)', 'actual G^{''}(s)', ...
'asympt. G^{"}(s)', 'actual comp. G^{"}(s)',...
'','actual no-comp. G^{"}(s)','Location', 'best');
legend(hax2, 'asympt. G^{''}(s)', 'actual G^{''}(s)', ...
'asympt. G^{"}(s)', 'actual comp. G^{"}(s)',...
'','actual no-comp. G^{"}(s)','Location', 'best');
Example 8 - Non Minimum Phase System: RHP Zero and Negative Static Gain
Consider the LTI systems described by the transfer functions:
Analyzing the transfer functions, you get

s = tf('s'); % the transfer function building helper
Gs1 = (1-s)/(s*(1+s))
Gs1 =
-s + 1
-------
s^2 + s
Continuous-time transfer function.
Model Properties
Gs2 = -(1-s)/(s*(1+s))
Gs2 =
s - 1
-------
s^2 + s
Continuous-time transfer function.
Model Properties
omVALS = logspace(-2,1, 1e4);
% We want to evaluate the frequency response at 10000
% angular pulsation values between 10^-2 and 10^1 rad/s.
% The greater the number of angular pulsation values
% used to evaluate the frequency response, the better
% the accuracy of the graphs of both the asymptotic
% let's create the diagrams: we need to store the handler of the figure
% refer to the help of drawBodediagrams( )
asBcolor = [1 0 0; 0 1 0]; % red or green color (RGB) for the asymptotic diagrams
Bcolors = [0 0.4470 0.7410; ...
0.9290 0.6940 0.1250; ...
0.4940 0.1840 0.5560; ...
0.4660 0.6740 0.1880; ...
0.3010 0.7450 0.9330; ...
0.8500 0.3250 0.0980; ...
0.6350 0.0780 0.1840]; % a different color for the actual diagram
hf = figure('Units','centimeters','Position',[0.01, 0.01, 24, 28]);
drawBodediagrams(Gs1, omVALS, asBcolor(1,:), 2.5, '-', ...
Bcolors(1,:), 1.5, '-', hf, true);
drawBodediagrams(Gs2, omVALS, asBcolor(2,:), 2.5, '-', ...
Bcolors(2,:), 2.5, '-', hf, true);
[hax1, hax2] = drawBodediagrams(Gs2, omVALS, asBcolor(2,:), 2.5, '-', ...
Bcolors(3,:), 1.5, '-', hf, false);
legend(hax1, 'asympt. G^{''}(s)', 'actual comp. G^{''}(s)', ...
'asympt. G^{"}(s)', 'actual comp. G^{"}(s)',...
'','actual no-comp. G^{"}(s)','Location', 'best');
legend(hax2, 'asympt. G^{''}(s)', 'actual comp. G^{''}(s)', ...
'asympt. G^{"}(s)', 'actual comp. G^{"}(s)',...
'','actual no-comp. G^{"}(s)','Location', 'best');
rmpath(genpath("BodeDiagram"))