From Lesson 1, the frequency response of an LTI system with transfer function \( G(s) \) is obtained by evaluation on the imaginary axis: \( G(j\omega) \), \( \omega \ge 0 \). A Bode plot is the pair of curves: \( 20\log_{10}|G(j\omega)| \) (magnitude in dB) and \( \angle G(j\omega) \) (phase in degrees), each vs. \( \log_{10}\omega \).

The central engineering idea is factorization into canonical building blocks (gains, integrators, first-order poles/zeros, second-order factors). Because logarithms turn products into sums, Bode plots can be constructed by adding the contributions of these factors.

flowchart TD
  A["Start with transfer function G(s)"] --> B["Factorize into standard blocks"]
  B --> C["Replace s by jw to get G(jw)"]
  C --> D["Magnitude: 20*log10(|G(jw)|)"]
  C --> E["Phase: arg(G(jw)) in degrees"]
  D --> F["Asymptotic straight-line sketch"]
  E --> G["Asymptotic phase sketch"]
  F --> H["Sum contributions (dB add)"]
  G --> I["Sum contributions (phase add)"]
  H --> J["Refine near corners if needed"]
  I --> J
        

Let \( G(j\omega) \in \mathbb{C} \). Write its polar form \( G(j\omega)=|G(j\omega)|e^{j\phi(\omega)} \) where \( \phi(\omega)=\arg(G(j\omega)) \). Define:

\[ M(\omega) \triangleq 20\log_{10}|G(j\omega)|,\quad \Phi(\omega) \triangleq \frac{180}{\pi}\arg(G(j\omega)). \]

The key construction identities follow from complex multiplication. Suppose \( G(s)=\prod_{k=1}^{N}G_k(s) \). Then for each fixed \( \omega \):

\[ |G(j\omega)| = \prod_{k=1}^{N}|G_k(j\omega)|, \qquad \arg(G(j\omega)) = \sum_{k=1}^{N}\arg(G_k(j\omega)) \; (\text{mod } 2\pi). \]

Proof (product magnitude and phase):

For complex numbers \( z_k = r_k e^{j\theta_k} \), the product is \( \prod_k z_k = \left(\prod_k r_k\right)e^{j\sum_k\theta_k} \). Hence the magnitude multiplies and the angles add (up to \( 2\pi \) wraps). Converting to Bode coordinates gives:

\[ 20\log_{10}|G(j\omega)| = 20\log_{10}\left(\prod_{k=1}^{N}|G_k(j\omega)|\right) = \sum_{k=1}^{N}20\log_{10}|G_k(j\omega)|. \]

Therefore, Bode magnitude (in dB) is additive across factors, and Bode phase is additive after consistent phase unwrapping.

Reality symmetry (important in practice). If \( G(s) \) has real coefficients, then \( G(-j\omega) = \overline{G(j\omega)} \), implying \( |G(-j\omega)|=|G(j\omega)| \) and \( \arg(G(-j\omega)) = -\arg(G(j\omega)) \) (up to wraps). This is why plots are typically drawn only for \( \omega \ge 0 \).

For hand construction, express the transfer function as a product of canonical factors. A typical real-rational form is:

\[ G(s) = K \frac{ \prod_{i=1}^{n_z}\left(1+\frac{s}{\omega_{z,i}}\right) \prod_{m=1}^{n_{z2}} \left(1 + 2\zeta_{z,m}\frac{s}{\omega_{n,z,m}} + \left(\frac{s}{\omega_{n,z,m}}\right)^2\right) }{ s^{n_0} \prod_{i=1}^{n_p}\left(1+\frac{s}{\omega_{p,i}}\right) \prod_{m=1}^{n_{p2}} \left(1 + 2\zeta_{p,m}\frac{s}{\omega_{n,p,m}} + \left(\frac{s}{\omega_{n,p,m}}\right)^2\right) }. \]

Here \( K \in \mathbb{R} \) is the gain, \( s^{n_0} \) represents \( n_0 \) integrators, \( \omega_{z,i},\omega_{p,i} \) are first-order break frequencies, and \( (\omega_{n},\zeta) \) parameterize second-order factors (introduced earlier in Chapter 10 as natural frequency and damping ratio).

flowchart LR
  G["G(s)"] --> K["Gain K"]
  G --> I["Integrators s^-n0"]
  G --> Z1["1st-order zeros (1+s/wz)"]
  G --> P1["1st-order poles 1/(1+s/wp)"]
  G --> Z2["2nd-order zeros (zeta, wn)"]
  G --> P2["2nd-order poles (zeta, wn)"]
        

We now derive the exact magnitude/phase for the most important factors used in construction rules.

(a) Constant gain \(K\)

\[ |K| = |K|,\quad M_K = 20\log_{10}|K|,\quad \Phi_K = \begin{cases} 0^\circ & \text{if } K > 0 \\ 180^\circ & \text{if } K < 0 \end{cases} \]

(b) Integrator \(1/s\) (and \(1/s^{n_0}\))

\[ \frac{1}{j\omega} = \frac{1}{\omega}e^{-j\pi/2} \;\Rightarrow\; \left|\frac{1}{j\omega}\right|=\frac{1}{\omega},\quad \arg\left(\frac{1}{j\omega}\right)=-\frac{\pi}{2}. \]

\[ M_{1/s}(\omega) = 20\log_{10}\left(\frac{1}{\omega}\right) = -20\log_{10}\omega,\quad \Phi_{1/s}(\omega) = -90^\circ. \]

For \( 1/s^{n_0} \), additivity yields: \( M(\omega) = -20n_0\log_{10}\omega \) and \( \Phi(\omega) = -90n_0^\circ \).

(c) First-order zero \( (1+s/\omega_c) \)

Substitute \( s=j\omega \) and let \( x=\omega/\omega_c \):

\[ 1 + jx = \sqrt{1+x^2}\;e^{j\arctan(x)}. \]

\[ \left|1+j\frac{\omega}{\omega_c}\right| = \sqrt{1+\left(\frac{\omega}{\omega_c}\right)^2}, \qquad \arg\left(1+j\frac{\omega}{\omega_c}\right)=\arctan\left(\frac{\omega}{\omega_c}\right). \]

(d) First-order pole \( 1/(1+s/\omega_c) \)

Since \( 1/(1+jx)= (1+jx)^{-1} \), magnitude inverts and phase negates:

\[ \left|\frac{1}{1+j\frac{\omega}{\omega_c}}\right| = \frac{1}{\sqrt{1+\left(\frac{\omega}{\omega_c}\right)^2}}, \qquad \arg\left(\frac{1}{1+j\frac{\omega}{\omega_c}}\right) = -\arctan\left(\frac{\omega}{\omega_c}\right). \]

Bode magnitude uses \( \log_{10}\omega \) on the horizontal axis. A key simplification is to approximate \( \sqrt{1+x^2} \) using its low/high-frequency limits.

First-order zero \( (1+s/\omega_c) \)

\[ M_z(\omega) = 20\log_{10}\sqrt{1+\left(\frac{\omega}{\omega_c}\right)^2} = 10\log_{10}\left(1+\left(\frac{\omega}{\omega_c}\right)^2\right). \]

For \( \omega/\omega_c \ll 1 \), \( 1+x^2 \approx 1 \): \( M_z(\omega)\approx 0 \) dB. For \( \omega/\omega_c \gg 1 \), \( 1+x^2 \approx x^2 \):

\[ M_z(\omega) \approx 10\log_{10}\left(\left(\frac{\omega}{\omega_c}\right)^2\right) = 20\log_{10}\left(\frac{\omega}{\omega_c}\right). \]

This yields the straight-line rule: 0 dB up to \( \omega_c \), then +20 dB/dec beyond \( \omega_c \).

Corner error (exact vs asymptote at \( \omega=\omega_c \))

\[ M_z(\omega_c) = 20\log_{10}\sqrt{2} = 10\log_{10}2 \approx 3.0103\;\text{dB}. \]

The asymptote gives \( 0 \) dB at \( \omega_c \), so the error is about \( +3.01 \) dB for a zero (and \( -3.01 \) dB for a pole).

First-order pole \( 1/(1+s/\omega_c) \)

\[ M_p(\omega) = -10\log_{10}\left(1+\left(\frac{\omega}{\omega_c}\right)^2\right) \approx \begin{cases} 0\;\text{dB} & \omega \ll \omega_c \\ -20\log_{10}\left(\frac{\omega}{\omega_c}\right)\;\text{dB} & \omega \gg \omega_c \end{cases} \]

Integrator \( 1/s \)

\[ M_{1/s}(\omega) = -20\log_{10}\omega \quad \Rightarrow \quad \text{slope } -20\;\text{dB/dec}. \]

Slope proof (general). If \( M(\omega)=20\log_{10}|G(j\omega)| \), define \( u=\log_{10}\omega \). Then:

\[ \frac{dM}{du} = 20\frac{d}{du}\log_{10}|G(j10^u)|. \]

For an integrator, \( |G|=\omega^{-1}=10^{-u} \), so \( M(u)=20\log_{10}(10^{-u})=-20u \), hence \( dM/du=-20 \) dB/dec exactly.

Time constant notation. Many texts write a first-order factor as \( (1 + s \tau) \) with \( \tau = 1/\omega_c \), so the break is \( \omega_c = 1/\tau \).

Phase is constructed by summing the phase contributions of each factor. For a first-order zero: \( \Phi_z(\omega)=\arctan(\omega/\omega_c)\cdot 180/\pi \). For a first-order pole: \( \Phi_p(\omega)=-\arctan(\omega/\omega_c)\cdot 180/\pi \).

Asymptotic phase approximation (standard “two-decade” rule)

For a first-order zero at \( \omega_c \), approximate phase as: \( 0^\circ \) up to \( 0.1\omega_c \), linearly ramp to \( +90^\circ \) at \( 10\omega_c \), then stay at \( +90^\circ \). The slope is \( +45^\circ/\text{dec} \) across that range. Poles are the negative of this.

Justification of the two-decade approximation

Evaluate the exact phase at the endpoints:

\[ \Phi_z(0.1\omega_c) = \frac{180}{\pi}\arctan(0.1) \approx 5.71^\circ, \qquad \Phi_z(10\omega_c) = \frac{180}{\pi}\arctan(10) \approx 84.29^\circ. \]

Thus “clamping” to \( 0^\circ \) and \( 90^\circ \) at those points introduces at most about \( 5.71^\circ \) endpoint error, while giving a rapid sketching rule.

Phase unwrapping

Because \( \arg(\cdot) \) is defined modulo \( 360^\circ \), computed phase may jump by \( \pm 360^\circ \). A continuous Bode phase uses unwrapping: add/subtract \( 360^\circ \) to remove discontinuities larger than \( 180^\circ \).

\[ \Phi_{\text{unwrap}}(\omega_k) \leftarrow \Phi(\omega_k) + 360^\circ \cdot m_k, \quad \text{choose integers } m_k \text{ to keep successive jumps } < 180^\circ. \]

Consider the transfer function (chosen to include a gain, an integrator, a first-order zero, and a first-order pole):

\[ G(s) = 10\;\frac{1 + \frac{s}{1}}{s\left(1+\frac{s}{10}\right)}. \]

The factors are: gain \( K=10 \), one integrator \( 1/s \), a zero at \( \omega_{z}=1 \) rad/s, and a pole at \( \omega_{p}=10 \) rad/s.

Asymptotic magnitude construction

• Gain: \( 20\log_{10}(10)=20 \) dB (vertical shift).
• Integrator: slope \( -20 \) dB/dec over all frequencies.
• Zero at \( 1 \): add \( +20 \) dB/dec for \( \omega \ge 1 \).
• Pole at \( 10 \): add \( -20 \) dB/dec for \( \omega \ge 10 \).

Therefore, slopes by region: \( \omega < 1 \): slope \( -20 \) dB/dec; \( 1 \le \omega < 10 \): slope \( 0 \) dB/dec; \( \omega \ge 10 \): slope \( -20 \) dB/dec.

Exact magnitude/phase expressions

\[ G(j\omega) = 10 \cdot \frac{1 + j\omega}{j\omega\left(1 + j\frac{\omega}{10}\right)}. \]

\[ |G(j\omega)| = 10\cdot \frac{\sqrt{1+\omega^2}}{\omega\sqrt{1+\left(\frac{\omega}{10}\right)^2}}, \qquad \Phi(\omega)= \arctan(\omega) - 90^\circ - \arctan\left(\frac{\omega}{10}\right). \]

These exact forms justify the additive slope changes and additive phase transitions used in the Bode rules.

This lab computes \( G(j\omega) \) directly, produces exact magnitude/phase curves, and overlays asymptotic approximations built from the rules in Sections 4–5. It also demonstrates library-based Bode plotting using python-control.

File: Chapter12_Lesson2.py

"""
Chapter 12 - Lesson 2: Bode Plots: Magnitude, Phase, Asymptotes, and Construction Rules
System Dynamics (Control Engineering)

This script demonstrates:
1) Library-based Bode plotting using python-control.
2) Manual computation of frequency response G(jw), magnitude (dB), phase (deg).
3) Construction of simple asymptotic Bode magnitude and phase approximations
   for first-order poles/zeros and integrators.

Dependencies (pip):
  - numpy
  - matplotlib
  - control  (python-control)
"""

from __future__ import annotations
import numpy as np
import matplotlib.pyplot as plt

try:
    import control  # type: ignore
except Exception as e:
    control = None
    print("python-control is not available:", e)


def polyval(coeffs: np.ndarray, s: np.ndarray) -> np.ndarray:
    """Evaluate polynomial with coefficients in descending powers."""
    y = np.zeros_like(s, dtype=np.complex128)
    for c in coeffs:
        y = y * s + c
    return y


def freq_response(num: np.ndarray, den: np.ndarray, w: np.ndarray) -> np.ndarray:
    """Compute G(jw) for a rational transfer function with real coefficients."""
    s = 1j * w
    return polyval(num, s) / polyval(den, s)


def mag_db(Gjw: np.ndarray) -> np.ndarray:
    return 20.0 * np.log10(np.abs(Gjw))


def unwrap_phase_deg(phase_deg: np.ndarray) -> np.ndarray:
    """Simple phase unwrap in degrees (adds/subtracts 360 to reduce jumps)."""
    out = phase_deg.copy()
    for k in range(1, len(out)):
        d = out[k] - out[k - 1]
        if d > 180.0:
            out[k:] -= 360.0
        elif d < -180.0:
            out[k:] += 360.0
    return out


def phase_deg(Gjw: np.ndarray) -> np.ndarray:
    ph = np.angle(Gjw, deg=True)
    return unwrap_phase_deg(ph)


def asym_mag_db_first_order_zero(w: np.ndarray, wc: float) -> np.ndarray:
    """Asymptotic magnitude (dB) for factor (1 + s/wc)."""
    y = np.zeros_like(w, dtype=float)
    mask = w >= wc
    y[mask] = 20.0 * np.log10(w[mask] / wc)
    return y


def asym_mag_db_first_order_pole(w: np.ndarray, wc: float) -> np.ndarray:
    """Asymptotic magnitude (dB) for factor 1/(1 + s/wc)."""
    return -asym_mag_db_first_order_zero(w, wc)


def asym_phase_deg_first_order_zero(w: np.ndarray, wc: float) -> np.ndarray:
    """
    Standard Bode phase approximation for (1 + s/wc):
      0 deg for w <= 0.1 wc
      linear ramp 0..+90 between 0.1 wc and 10 wc
      +90 deg for w >= 10 wc
    """
    y = np.zeros_like(w, dtype=float)
    w1, w2 = 0.1 * wc, 10.0 * wc
    mid = (w > w1) & (w < w2)
    y[w >= w2] = 90.0
    y[mid] = 45.0 * (np.log10(w[mid] / w1))
    return y


def asym_phase_deg_first_order_pole(w: np.ndarray, wc: float) -> np.ndarray:
    return -asym_phase_deg_first_order_zero(w, wc)


def asym_mag_db_integrator(w: np.ndarray, m: int = 1) -> np.ndarray:
    """Asymptotic magnitude for 1/s^m: -20m dB/dec relative to w=1 rad/s."""
    return -20.0 * m * np.log10(w)


def asym_phase_deg_integrator(m: int = 1) -> float:
    return -90.0 * m


# Example transfer function:
#   G(s) = 10 * (1 + s/1) / ( s * (1 + s/10) )
num = np.array([10.0, 10.0])
den = np.array([0.1, 1.0, 0.0])

w = np.logspace(-2, 3, 2000)  # rad/s
Gjw = freq_response(num, den, w)
M = mag_db(Gjw)
P = phase_deg(Gjw)

# Asymptotes: gain 10 (20 dB), one zero at wc=1, one pole at wc=10, one integrator
Mag_asym = 20.0 * np.log10(10.0) \
    + asym_mag_db_first_order_zero(w, wc=1.0) \
    + asym_mag_db_first_order_pole(w, wc=10.0) \
    + asym_mag_db_integrator(w, m=1)

Phase_asym = 0.0 \
    + asym_phase_deg_first_order_zero(w, wc=1.0) \
    + asym_phase_deg_first_order_pole(w, wc=10.0) \
    + asym_phase_deg_integrator(m=1)

plt.figure()
plt.semilogx(w, M, label="Exact magnitude (dB)")
plt.semilogx(w, Mag_asym, "--", label="Asymptotic magnitude (dB)")
plt.xlabel("Frequency w [rad/s]"); plt.ylabel("Magnitude [dB]")
plt.grid(True, which="both"); plt.legend()

plt.figure()
plt.semilogx(w, P, label="Exact phase (deg)")
plt.semilogx(w, Phase_asym, "--", label="Asymptotic phase (deg)")
plt.xlabel("Frequency w [rad/s]"); plt.ylabel("Phase [deg]")
plt.grid(True, which="both"); plt.legend()

if control is not None:
    sys_tf = control.TransferFunction(num, den)
    plt.figure()
    control.bode_plot(sys_tf, w, dB=True, Hz=False, deg=True, plot=True)

plt.show()
      

Related libraries for System Dynamics work in Python include: python-control (transfer functions, frequency response), scipy.signal (LTI utilities), and numpy (complex evaluation).

MATLAB’s Control System Toolbox provides exact Bode plots via tf and bode, while the underlying definitions can be reproduced by direct evaluation of \( G(j\omega) \).

File: Chapter12_Lesson2.m

% Chapter 12 - Lesson 2: Bode Plots: Magnitude, Phase, Asymptotes, and Construction Rules
% System Dynamics (Control Engineering)

clear; clc; close all;

% Example transfer function:
%   G(s) = 10 * (1 + s/1) / ( s * (1 + s/10) )
K  = 10;
wz = 1;    % zero break (rad/s)
wp = 10;   % pole break (rad/s)

s = tf('s');
G = K * (1 + s/wz) / ( s * (1 + s/wp) );

w = logspace(-2, 3, 2000);  % rad/s

% Exact Bode using toolbox
figure;
bode(G, w); grid on;
title('Exact Bode (MATLAB bode)');

% Manual frequency response
jw = 1i*w;
Gj = K * (1 + jw/wz) ./ ( jw .* (1 + jw/wp) );

Mag = 20*log10(abs(Gj));
Ph  = unwrap(angle(Gj)) * 180/pi;

% Asymptotic magnitude
MagAsym = 20*log10(K) ...
    + (w >= wz) .* (20*log10(w/wz)) ...
    - (w >= wp) .* (20*log10(w/wp)) ...
    - 20*log10(w);  % integrator 1/s

% Asymptotic phase (two-decade rule)
phaseZero = zeros(size(w));
phasePole = zeros(size(w));

w1 = 0.1*wz; w2 = 10*wz;
mid = (w > w1) & (w < w2);
phaseZero(w >= w2) = 90;
phaseZero(mid) = 45*log10(w(mid)/w1);

w1 = 0.1*wp; w2 = 10*wp;
mid = (w > w1) & (w < w2);
phasePole(w >= w2) = -90;
phasePole(mid) = -45*log10(w(mid)/w1);

PhAsym = phaseZero + phasePole - 90;

figure;
semilogx(w, Mag, 'LineWidth', 1.2); hold on;
semilogx(w, MagAsym, '--', 'LineWidth', 1.2);
grid on; xlabel('w [rad/s]'); ylabel('Magnitude [dB]');
legend('Exact', 'Asymptotic'); title('Magnitude: exact vs asymptotic');

figure;
semilogx(w, Ph, 'LineWidth', 1.2); hold on;
semilogx(w, PhAsym, '--', 'LineWidth', 1.2);
grid on; xlabel('w [rad/s]'); ylabel('Phase [deg]');
legend('Exact', 'Asymptotic'); title('Phase: exact vs asymptotic');
      

For a model-based workflow, the following companion script builds a Simulink model of the transfer function and attempts linearization for Bode plotting.

File: Chapter12_Lesson2_Simulink.m

% Chapter 12 - Lesson 2 (Simulink companion): Bode plot from a Simulink model

clear; clc;

model = 'Chapter12_Lesson2_BodeModel';
if bdIsLoaded(model)
    close_system(model, 0);
end
new_system(model);
open_system(model);

add_block('simulink/Sources/In1', [model '/u'], 'Position', [30 80 60 100]);
add_block('simulink/Continuous/Transfer Fcn', [model '/G'], 'Position', [120 70 230 110]);
add_block('simulink/Sinks/Out1', [model '/y'], 'Position', [300 80 330 100]);

% Example G(s) = 10*(s+1) / (0.1*s^2 + s)
set_param([model '/G'], 'Numerator', '[10 10]', 'Denominator', '[0.1 1 0]');

add_line(model, 'u/1', 'G/1');
add_line(model, 'G/1', 'y/1');

set_param(model, 'StopTime', '10');

try
    io(1) = linio([model '/u'], 1, 'input');
    io(2) = linio([model '/y'], 1, 'output');
    sys = linearize(model, io);
    figure; bode(sys); grid on;
    title('Bode from linearized Simulink model');
catch ME
    warning('Linearization could not be performed. Details:\n%s', ME.message);
    disp('Alternative: Use Simulink "Linear Analysis" tool to define I/O and plot Bode.');
end
      

In embedded and real-time contexts, you may compute frequency response numerically (complex polynomial evaluation) and export results for plotting. The implementations below avoid control toolboxes and implement the definition directly.

File: Chapter12_Lesson2.cpp

/*
Chapter 12 - Lesson 2: Bode Plots: Magnitude, Phase, Asymptotes, and Construction Rules
System Dynamics (Control Engineering)
*/

#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
#include <fstream>
#include <iomanip>

static std::complex<double> polyval(const std::vector<double>& c,
                                    const std::complex<double>& s) {
    std::complex<double> y(0.0, 0.0);
    for (double a : c) {
        y = y * s + a;
    }
    return y;
}

static double mag_db(const std::complex<double>& z) {
    return 20.0 * std::log10(std::abs(z));
}

static double phase_deg(const std::complex<double>& z) {
    return std::atan2(z.imag(), z.real()) * 180.0 / M_PI;
}

static void unwrap_phase(std::vector<double>& ph_deg) {
    for (size_t k = 1; k < ph_deg.size(); ++k) {
        double d = ph_deg[k] - ph_deg[k - 1];
        if (d > 180.0) {
            for (size_t j = k; j < ph_deg.size(); ++j) ph_deg[j] -= 360.0;
        } else if (d < -180.0) {
            for (size_t j = k; j < ph_deg.size(); ++j) ph_deg[j] += 360.0;
        }
    }
}

int main() {
    // Example: G(s) = 10*(1 + s/1) / ( s*(1 + s/10) )
    const std::vector<double> num = {10.0, 10.0};
    const std::vector<double> den = {0.1, 1.0, 0.0};

    const int N = 2000;
    const double w_min = 1e-2;
    const double w_max = 1e3;

    std::vector<double> w(N), mag(N), ph(N);

    for (int k = 0; k < N; ++k) {
        double t = static_cast<double>(k) / static_cast<double>(N - 1);
        double wk = w_min * std::pow(w_max / w_min, t); // logspace
        w[k] = wk;

        std::complex<double> s(0.0, wk); // jw
        std::complex<double> Gjw = polyval(num, s) / polyval(den, s);

        mag[k] = mag_db(Gjw);
        ph[k]  = phase_deg(Gjw);
    }

    unwrap_phase(ph);

    std::ofstream f("bode_output.csv");
    f << "w,mag_db,phase_deg\n";
    f << std::setprecision(12);
    for (int k = 0; k < N; ++k) {
        f << w[k] << "," << mag[k] << "," << ph[k] << "\n";
    }
    f.close();

    std::cout << "Wrote bode_output.csv (" << N << " points)." << std::endl;
    return 0;
}
      

File: Chapter12_Lesson2.java

/*
Chapter 12 - Lesson 2: Bode Plots: Magnitude, Phase, Asymptotes, and Construction Rules
System Dynamics (Control Engineering)
*/

import java.io.PrintWriter;
import java.util.Locale;

public class Chapter12_Lesson2 {

    static final class Complex {
        final double re;
        final double im;

        Complex(double re, double im) { this.re = re; this.im = im; }

        Complex add(Complex other) {
            return new Complex(this.re + other.re, this.im + other.im);
        }

        Complex mul(Complex other) {
            return new Complex(this.re * other.re - this.im * other.im,
                               this.re * other.im + this.im * other.re);
        }

        Complex div(Complex other) {
            double denom = other.re * other.re + other.im * other.im;
            return new Complex((this.re * other.re + this.im * other.im) / denom,
                               (this.im * other.re - this.re * other.im) / denom);
        }

        double abs() { return Math.hypot(this.re, this.im); }

        double phaseDeg() { return Math.atan2(this.im, this.re) * 180.0 / Math.PI; }
    }

    static Complex polyval(double[] c, Complex s) {
        Complex y = new Complex(0.0, 0.0);
        for (double a : c) {
            y = y.mul(s).add(new Complex(a, 0.0));
        }
        return y;
    }

    static double magDb(Complex z) { return 20.0 * Math.log10(z.abs()); }

    static void unwrapPhase(double[] phDeg) {
        for (int k = 1; k < phDeg.length; ++k) {
            double d = phDeg[k] - phDeg[k - 1];
            if (d > 180.0) {
                for (int j = k; j < phDeg.length; ++j) phDeg[j] -= 360.0;
            } else if (d < -180.0) {
                for (int j = k; j < phDeg.length; ++j) phDeg[j] += 360.0;
            }
        }
    }

    public static void main(String[] args) throws Exception {
        Locale.setDefault(Locale.ROOT);

        double[] num = new double[]{10.0, 10.0};
        double[] den = new double[]{0.1, 1.0, 0.0};

        int N = 2000;
        double wMin = 1e-2;
        double wMax = 1e3;

        double[] w = new double[N];
        double[] mag = new double[N];
        double[] ph = new double[N];

        for (int k = 0; k < N; ++k) {
            double t = (double) k / (double) (N - 1);
            double wk = wMin * Math.pow(wMax / wMin, t); // logspace
            w[k] = wk;

            Complex s = new Complex(0.0, wk); // jw
            Complex Gjw = polyval(num, s).div(polyval(den, s));

            mag[k] = magDb(Gjw);
            ph[k] = Gjw.phaseDeg();
        }

        unwrapPhase(ph);

        try (PrintWriter out = new PrintWriter("bode_output_java.csv")) {
            out.println("w,mag_db,phase_deg");
            for (int k = 0; k < N; ++k) {
                out.printf(Locale.ROOT, "%.12g,%.12g,%.12g%n", w[k], mag[k], ph[k]);
            }
        }
    }
}
      

Typical C++/Java “system dynamics” support libraries are less standardized than Python/MATLAB; therefore, direct evaluation via complex arithmetic is an important baseline method. When available, common numerical back-ends include Eigen (C++) and Apache Commons Math (Java) for complex/vector utilities.

Mathematica natively supports symbolic and numeric transfer-function models. The notebook below defines the same example and produces a Bode plot, then computes \( G(j\omega) \) directly.

File: Chapter12_Lesson2.nb

(* Chapter 12 - Lesson 2: Bode Plots (Wolfram Mathematica) *)

(* Define transfer function: G(s)=10 (1+s)/(s (1+s/10)) *)
G = TransferFunctionModel[10 (1 + s)/(s (1 + s/10)), s];

w = 10^Range[-2, 3, 0.0025];

BodePlot[G, {w[[1]], w[[-1]]}, ScalingFunctions -> "Log", PlotLegends -> None]

(* Manual frequency response and asymptotes *)
jw = I w;
Gj = 10 (1 + jw)/(jw (1 + jw/10));
Mag = 20 Log10[Abs[Gj]];
Ph = Arg[Gj];

ListLinePlot[Transpose[{w, Mag}], ScalingFunctions -> {"Log", None}, AxesLabel -> {"w", "Mag (dB)"}]
      

Problem 1 (Corner error for first-order factors): Let \( H(s)=1+\frac{s}{\omega_c} \). Show that the asymptotic magnitude at \( \omega=\omega_c \) differs from the exact magnitude by \( 20\log_{10}\sqrt{2} \) dB, and compute the numerical value.

Solution: Evaluate exactly at \( s=j\omega_c \):

\[ |H(j\omega_c)| = \left|1 + j\right| = \sqrt{2} \;\Rightarrow\; 20\log_{10}|H(j\omega_c)| = 20\log_{10}\sqrt{2} = 10\log_{10}2 \approx 3.0103\;\text{dB}. \]

The standard magnitude asymptote for a zero is \( 0 \) dB at the corner, hence the asymptote underestimates by \( 3.0103 \) dB. For a pole, the sign reverses.

Problem 2 (Hand-sketch slopes): Construct the asymptotic magnitude slope (in dB/dec) of \( G(s)=\frac{K}{s(1+s/\omega_p)} \) for \( \omega < \omega_p \) and \( \omega \ge \omega_p \). Do not compute intercepts; report slopes only.

Solution:

Factor contributions: \( K \) adds only a constant shift (0 slope), \( 1/s \) contributes \( -20 \) dB/dec at all frequencies, and \( 1/(1+s/\omega_p) \) contributes \( 0 \) dB/dec for \( \omega < \omega_p \) and \( -20 \) dB/dec for \( \omega \ge \omega_p \).

Hence: \( \omega < \omega_p \Rightarrow \) slope \( -20 \) dB/dec, \( \omega \ge \omega_p \Rightarrow \) slope \( -40 \) dB/dec.

Problem 3 (Phase additivity and integrator contribution): For \( G(s)=10\frac{1+s}{s(1+s/10)} \), show that \( \Phi(\omega) = \arctan(\omega) - 90^\circ - \arctan(\omega/10) \).

Solution:

Substitute \( s=j\omega \): \( G(j\omega)=10\cdot (1+j\omega)/\big(j\omega(1+j\omega/10)\big) \). Using \( \arg(ab)=\arg(a)+\arg(b) \) and \( \arg(1/a)=-\arg(a) \):

\[ \arg(G(j\omega)) = \arg(10) + \arg(1+j\omega) - \arg(j\omega) - \arg\left(1+j\frac{\omega}{10}\right). \]

For \( \omega > 0 \), \( \arg(10)=0^\circ \), \( \arg(1+j\omega)=\arctan(\omega)\cdot 180/\pi \), \( \arg(j\omega)=90^\circ \), and \( \arg(1+j\omega/10)=\arctan(\omega/10)\cdot 180/\pi \). Therefore:

\[ \Phi(\omega) = \arctan(\omega) - 90^\circ - \arctan\left(\frac{\omega}{10}\right). \]

Problem 4 (Reality symmetry): Assume \( G(s)=\frac{N(s)}{D(s)} \) where \( N,D \) have real coefficients. Prove that \( G(-j\omega)=\overline{G(j\omega)} \) for real \( \omega \), and conclude that magnitude is even and phase is odd (up to wrapping).

Solution:

For any polynomial with real coefficients, \( \overline{N(j\omega)} = N(\overline{j\omega}) = N(-j\omega) \). Similarly, \( \overline{D(j\omega)} = D(-j\omega) \). Thus:

\[ \overline{G(j\omega)} = \overline{\frac{N(j\omega)}{D(j\omega)}} = \frac{\overline{N(j\omega)}}{\overline{D(j\omega)}} = \frac{N(-j\omega)}{D(-j\omega)} = G(-j\omega). \]

Taking magnitudes gives \( |G(-j\omega)| = |\overline{G(j\omega)}|=|G(j\omega)| \). Taking arguments yields \( \arg(G(-j\omega))=\arg(\overline{G(j\omega)})=-\arg(G(j\omega)) \) (mod \( 2\pi \)), which explains the even/odd symmetry commonly used when plotting only \( \omega \ge 0 \).

Bode plots encode the frequency response \( G(j\omega) \) as logarithmic magnitude (dB) and unwrapped phase (deg). By factorizing \( G(s) \) into canonical blocks, we proved the additivity rules that enable rapid construction. We derived exact magnitude/phase formulas for integrators and first-order poles/zeros, developed straight-line asymptotes, and quantified corner errors (notably the \( \pm 3.01 \) dB corner deviation for first-order factors). These tools are prerequisites for resonance and bandwidth analysis in the next lesson.