Frequency Response & Bode Diagrams

Table of Contents

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
block_diagramGs.png
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

block_diagramGs.png
Then
What about using MATLAB?
clear
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)
magG2 = 0.4472
disp(sqrt(5)/5)
0.4472
% the phase angle of FreqResp_Gs_2 <--> arg G(j2) <-- NB rad
argG2_rad = angle(FreqResp_Gs_2)
argG2_rad = -1.1071
% converting from rad to deg
argG2_deg = rad2deg(argG2_rad)
argG2_deg = -63.4349
atan(2)
ans = 1.1071

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:
clear
 
xiVAL = [0.05, 0.2, 0.5, 0.7, 0.9]; % the damping factor values
Nxi = numel(xiVAL); % how many damping factors?
 
Gs_xi = cell(Nxi, 1);
% 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)
for n=1:Nxi
denGs = [ (1) (20*xiVAL(n)) 100];
% the coeffs of the polynomial at the denominator in G(s),
% in descending order
Gs_xi{n} = tf(numGs, denGs); % let's define the transfer function G(s)
% and store it in the data-structure
end % for
 
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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]);
 
for n=1:Nxi
[hax1, hax2] = drawBodediagrams(Gs_xi{n}, omVALS, asBcolor, 2.5, '-', ...
Bcolors(n, :), 1.5, '-', hf);
end % for
legend(hax1, 'asympt.', ['$\xi = ',num2str(xiVAL(1))], ...
'', ['\xi = ',num2str(xiVAL(2))], ...
'', ['\xi = ',num2str(xiVAL(3))], ...
'', ['\xi = ',num2str(xiVAL(4))], ...
'', ['\xi = ',num2str(xiVAL(5))], ...
'Location', 'best');
 
legend(hax2, 'asympt.', ['\xi = ',num2str(xiVAL(1))], ...
'', ['\xi = ',num2str(xiVAL(2))], ...
'', ['\xi = ',num2str(xiVAL(3))], ...
'', ['\xi = ',num2str(xiVAL(4))], ...
'', ['\xi = ',num2str(xiVAL(5))], ...
'Location', 'best');
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
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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, '-', ...
Bcolor, 1.5, '-', hf);
 
legend(hax1, 'asympt.', 'actual', ...
'Location', 'best');
 
legend(hax2, 'asympt.', 'actual', ...
'Location', 'best');
 

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
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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, '-', ...
% Bcolor, 1.5, '-', hf);
%
% legend(hax1, 'asympt.', 'actual', ...
% 'Location', 'best');
%
% legend(hax2, 'asympt.', 'actual', ...
% 'Location', 'best');

Example 4: an Unstable RHP Pole

Consider the LTI system described by the transfer function:
Analyzing the transfer function, you get
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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, '-', ...
Bcolor, 1.5, '-', hf);
 
legend(hax1, 'asympt.', 'actual', ...
'Location', 'best');
 
legend(hax2, 'asympt.', 'actual', ...
'Location', 'best');

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
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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)',...
'Location', 'best');
 
legend(hax2, 'asympt.', 'actual G^{''}(s)',...
'', 'actual G^{"}(s)', ...
'Location', 'best');
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
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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
clear
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.
%
% ** Note **:
% 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
% and actual diagrams.
 
% 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"))