Dynamical Systems:

Forced Response

Table of Contents





1. Linear Time-Invariant (LTI) Systems

The system $H$ represents a transformation (or rule) that maps an input signal $x(t)$ to an output signal $y(t)$:


$$ y(t) = H\{x(t)\}. $$



System Examples:


$$ \begin{aligned} \text{Identity:} &\quad y(t) = x(t) &\quad &\forall t \\ \text{Scaling:} &\quad y(t) = 2x(t) &\quad &\forall t \\ \text{Offset:} &\quad y(t) = x(t) + 2 &\quad &\forall t \\ \text{Square Signal:} &\quad y(t) = x(t)^2 &\quad &\forall t \\ \text{Shift:} &\quad y(t) = x(t+2) &\quad &\forall t \\ \text{Decimate:} &\quad y(t) = x(2t) &\quad &\forall t \\ \text{Square Time:} &\quad y(t) = x(t^2) &\quad &\forall t \end{aligned} $$


1.1. Linear Systems

A system $H$ is called linear if it satisfies:


(1) Scaling:


$$ H\{\alpha x(t)\} = \alpha H\{x(t)\}, \quad \forall \alpha \in \mathbb{C} $$



(2) Additivity:


If $y_1 = H{x_1}$ and $y_2 = H{x_2}$, then:


$$ H\{x_1 + x_2\} = y_1 + y_2. $$


1.2. Time-Invariant Systems

A system $H$ is time-invariant if shifting the input signal causes an identical shift in the output:


$$ \text{If } y(t) = H\{x(t)\}, \quad \text{then } H\{x(t - \tau)\} = y(t - \tau), \quad \forall \tau. $$



1.3. Linear Time-Invariant (LTI) Systems

We focus on systems that are both linear and time-invariant:


$$ H\{\alpha_1 x_1 + \alpha_2 x_2\}(t) = \alpha_1 H\{x_1\}(t) + \alpha_2 H\{x_2\}(t), $$


and:


$$ H\{x(t-\tau)\} = H\{x\}(t-\tau). $$


$$\begin{align*} \text{Identity} &\qquad y(t) = x(t) \quad &\forall t &\qquad \text{Linear} &\qquad \text{Time Invariant}\\ \text{Scaling} &\qquad y(t) = 2\,x(t) \quad &\forall t &\qquad \text{Linear} &\qquad \text{Time Invariant}\\ \text{Offset} &\qquad y(t) = x(t)\, + \,2 \quad &\forall t &\qquad \text{Non Linear} &\qquad \text{Time Invariant}\\ \text{Square signal} &\qquad y(t) = (x(t))^2 \quad &\forall t &\qquad \text{Non Linear} &\qquad \text{Time Invariant}\\ \text{Shift} &\qquad y(t) = x(t + 2) \quad &\forall t &\qquad \text{Linear} &\qquad \text{Time Invariant}\\ \text{Decimate} &\qquad y(t) = x(2t) \quad &\forall t &\qquad \text{Linear} &\qquad \text{Time Variant}\\ \text{Square time} &\qquad y(t) = x(t^2) \quad &\forall t &\qquad \text{Linear} &\qquad \text{Time Variant}\\ \end{align*}$$





2. Response to Non-Zero Input


So far, we have examined the natural response of systems with zero input and non-zero initial conditions:


$$ \begin{aligned} \dot{y} + \frac{1}{\tau}y &= 0 \\\\ \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2y &= 0 \end{aligned} $$



Response to non-zero constant input

Now, let’s analyze the response to non-zero constant input, which leads to inhomogeneous differential equations.



Assume the system is stable.


For constant forcing terms $q'$ or $p'$, the equations become:


$$ \begin{aligned} \dot{y} + \frac{1}{\tau}y &= q' \\\\ \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2y &= p' \end{aligned} $$


The system retains the same dynamics but converges to a different steady state.

Rewriting the equations by defining shifted variables:


$$ \begin{aligned} (\dot{y}-\dot{q}) + \frac{1}{\tau}(y-q) &= 0 \\\\ (\ddot{y}-\ddot{p}) + 2\zeta\omega_n(\dot{y}-\dot{p}) + \omega_n^2(y-p) &= 0 \end{aligned} $$


This shows that:

  • Transient response appears immediately after system activation and decays if the system is stable.
  • Steady-state response dominates in the long run, with all transients vanishing.

In steady state:


$$ \frac{dy}{dt} \to 0, \quad \frac{d^2y}{dt^2} \to 0 $$



Example 1: First-Order System


$$ \dot{y} + \frac{1}{\tau}y = q, \quad y(0) = 0 $$


At steady state:


$$ \dot{y}(\infty) = 0 \implies y(\infty) = \tau q $$



Example 2: Second-Order System


$$ \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2y = p, \quad y(0) = 0, \quad \dot{y}(0) = 0 $$


At steady state:


$$ \ddot{y}(\infty) = \dot{y}(\infty) = 0 \implies y(\infty) = \frac{p}{\omega_n^2} $$



Example 3: Mass-Spring-Damper System


(1) Horizontal Setup:


$$ \ddot{y} + \frac{c}{m}\dot{y} + \frac{k}{m}y = 0, \quad y(0) = y_0, \quad \dot{y}(0) = 0 $$



(2) Vertical Setup (subject to gravity):

  • $y(0) = 0 \;$ no initial displacement
  • $\dot{y}(0) = 0\;$ initially at rest

$$ \begin{aligned} m\ddot{y} + c\dot{y} + ky &= mg \\\\ \ddot{y} + \frac{c}{m}\dot{y} + \frac{k}{m}y &= g, \quad y(0) = 0, \quad \dot{y}(0) = 0 \end{aligned} $$



By shifting the origin to the static equilibrium point:


$$ y_{\text{eq}} = \frac{mg}{k} $$


the system behaves as a natural response starting from:


$$ y(0) = -\frac{mg}{k}, \quad \dot{y}(0) = 0 $$


General Interpretation:


$$ \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2y = \omega_n^2 p $$


Rewriting:


$$ (\ddot{y}-\ddot{p}) + 2\zeta\omega_n(\dot{y}-\dot{p}) + \omega_n^2(y-p) = 0 $$


This shows the system is a shifted version of the natural response around $p$.





3. Time Response to General Inputs

Up to this point, we have examined system responses to constant inputs or initial conditions. Now we move to general inputs and how an LTI system behaves in these scenarios. The main tool we use is the convolution integral.


3.1. Step Response: Example and Interpretation

Consider a simple first-order system:


$$ \dot{x} + 5x = 1 \quad \text{for} \quad t \geq 0, \qquad x(0) = 0$$


$$ \dot{x} + 5x = u(t), \quad x(0)=0, $$


where $u(t)$ is the unit step function:


$$ u(t) = \begin{cases} 1, & t \geq 0, \\ 0, & t < 0. \end{cases} $$



This models a system that is initially at rest, and at time $t=0$, a constant input of 1 is suddenly applied. The solution is:


$$ x(t) = \frac{1}{5}(1 - e^{-5t}). $$


  • The system starts at zero and asymptotically approaches the steady-state value $x(\infty)=\frac{1}{5}$.
  • The transient term $e^{-5t}$ decays exponentially with a time constant $\tau=1/5$.

In [ ]:
t = linspace(0,2,100);
x = 1/5*(1-exp(-5*t));

plot(t,x,t,0.2*ones(size(t)),'k--')
ylim([0,0.25])
xlabel('time')

No description has been provided for this image

3.2. Impulse Response

The Unit Impulse Signal: Dirac Delta Function

The Dirac delta function $\delta(t)$ serves as a mathematical model of an idealized pulse—zero width, infinite height, and unit area. This is a conceptual construct used to model instantaneous impulses applied at a specific point in time, such as a sudden jolt of force or current.


Formal Definition


$$\delta(t) = \lim_{\epsilon\rightarrow0}p_\epsilon(t)$$


where $p_\epsilon(t)$ is a rectangular pulse of width $2\epsilon$ and height $\frac{1}{2\epsilon}$, having a unit area regardless of $\epsilon$:


$$ \int_{-\epsilon}^{\epsilon} p_\epsilon(t) \, dt = 1. $$



As $\epsilon \to 0$, the pulse becomes an infinitely thin spike at $t=0$ with unit area—this is $\delta(t)$.

The unit-impulse function is represented by an arrow with the number 1, which represents its area.



Key Properties of $\delta(t)$


(1) Zero everywhere except at $t=0$:


$$ \delta(t) = \begin{cases} \infty, & t=0, \\[8pt] 0, & t \ne 0. \end{cases} $$


(2) Unit area under the curve:


$$ \int_{-\infty}^{\infty} \delta(t) \, dt = 1. $$


(3) Sifting property: It "samples" functions at a point:


$$ \begin{align*} \int_{-\infty}^{\infty}x(t)\delta(t)dt &= \int_{-\infty}^{0^-}x(t)\delta(t)dt + \int_{0^-}^{0^+}x(t)\delta(t)dt + \int_{0^+}^{\infty}x(t)\delta(t)dt \\&= 0 + x(0)\int_{0^-}^{0^+}\delta(t)dt + 0 \\&= x(0) \end{align*}$$


or more generally:


$$ \int_{-\infty}^{\infty} x(\tau) \delta(t-\tau) \, d\tau = x(t). $$



3.2.1. Impulse and Initial Velocity in Mass-Spring-Damper Systems


1. Question: How to Realize an Initial Velocity $v_0 \ne 0$?



The key idea comes from Newton's Second Law and the concept of impulse:


  • Impulse is the integral of force over time:

$$ \Delta p = \int_{0^-}^{0^+} F(t) dt = m v_0 - m(0) = m v_0. $$


  • A delta function force $F(t) = m v_0 \delta(t)$ applied at $t=0$ instantaneously changes the momentum of the mass:

$$ F(t) = m v_0 \delta(t), $$


$\quad$ which results in:


$$ \dot{y}(0^+) = v_0, \quad y(0) = 0. $$


This models an instantaneous "kick" to the mass.


2. Mathematical Model

For a mass-spring-damper system:


$$ m\ddot{y} + c\dot{y} + k y = F(t). $$


With an impulse force:


$$ F(t) = m v_0 \delta(t), $$


this becomes:


$$ m\ddot{y} + c\dot{y} + k y = m v_0 \delta(t), $$


or after dividing by $m$:


$$ \ddot{y} + \frac{c}{m}\dot{y} + \frac{k}{m}y = v_0 \delta(t). $$


The impulse at $t=0$ injects momentum:


$$ \Delta p = m v_0. $$


This results in:


$$ v_0 \delta(t) \quad \longrightarrow \quad \begin{cases} y(0) = 0,\\[10pt] \dot{y}(0) = v_0, \end{cases} \quad \longrightarrow \quad \text{natural response}. $$


Interpretation:

After the impulse (a “kick”), the system behaves as if it started from rest at $y=0$ but with an initial velocity $\dot{y}(0)=v_0$. This natural response evolves based on the system's inherent dynamics.



3.2.2. Impulse Response to LTI System



The impulse response is a fundamental concept for analyzing Linear Time-Invariant (LTI) systems. It characterizes how a system reacts to a very short input of unit area, namely the Dirac delta function $\delta(t)$.


Consider the second-order system:


$$ \ddot{y} + 2\zeta \omega_n \dot{y} + \omega_n^2 y = \delta(t) $$


where:

$$ \delta(t): \text{input}, \quad y(t): \text{output} $$


The response $y(t)$ to $\delta(t)$ fully characterizes the system's behavior. The impulse response $h(t)$ can later be used to compute the response to any input $x(t)$ via convolution:


$$ y(t) = h(t) * x(t) = \int_{-\infty}^{\infty} h(\tau) x(t-\tau) d\tau $$


Example:

Impulse Response of a First-Order System


$$ \dot{y} + 5y = \delta(t), \quad y(0) = 0 $$


The solution to this equation is:


$$ y(t) = h(t) = e^{-5t}, \quad t \ge 0 $$


Why is this the solution?

  • Interpretation: The system experiences an impulse $\delta(t)$ at $t=0$, causing an instantaneous change in $y$.
  • Impulse integral:

$$ \int_{0^-}^{0^+} \delta(t) dt = u(0^+) - u(0^-) = 1 $$


which can be interpreted as an instantaneous jump of 1 in the output.


After the impulse at $t=0$, the system evolves according to its homogeneous equation:


$$ \dot{y} + 5y = 0, \quad y(0) = 1 $$


leading to the natural response:


$$ y(t) = e^{-5t} $$


Key Insight

  • The impulse response is the system's fingerprint.

  • It reveals how the system reacts to the simplest possible input.

  • It forms the foundation for determining the system's output to any input via convolution.


In [ ]:
t = linspace(0,2,100);
h = exp(-5*t);

plot(t,h,t,zeros(size(t)), 'k--'),
ylim([-0.1,1.1])
xlabel('time')

No description has been provided for this image

3.3. Step Response

Consider the system:


$$ \begin{align*} \ddot{x} + 2\zeta\omega_n\dot{x}+\omega^2_nx &= u(t) &\text{or} \\ \\ \ddot{x} + 2\zeta\omega_n\dot{x}+\omega^2_nx &= 1, &t > 0 \end{align*} $$


The solution has two components:

  • Transient response, determined by natural modes (eigenvalues of the system).
  • Steady-state response, driven by the input $u(t)$.


Key point: The derivative of the step response $y_2(t)$ gives the impulse response $y_1(t)$.



3.4. Response to a General Input (in Time)

So far, we've analyzed systems with constant or simple inputs (such as step and impulse). But real-world systems often encounter more complex inputs $x(t)$. How can we determine the system’s response $y(t)$ for any input?


Setup:

Consider the first-order linear differential equation:


$$ \dot{y} + 5 y = x(t), \quad y(0) = 0 $$


This represents an LTI system (linear, time-invariant).


Key Idea:

The solution can be expressed using the system's impulse response $h(t)$:


$$ y(t) = h(t) * x(t) $$


where the symbol $*$ denotes convolution.



If this is true, we can compute output response to any general input if an impulse response is given


Convolution Definition:

For continuous-time signals:


$$ \begin{align*} y(t) = (h * x)(t) &= \int_{-\infty}^{\infty} h(\tau) x(t - \tau) d\tau \\ &= \int_{-\infty}^{\infty} x(\tau) h(t - \tau) d\tau \end{align*} $$


Why This Works:

Time-Invariance and Linearity


(1) Time-invariance ensures that if an impulse input $\delta(t-\tau)$ is shifted by $\tau$, the system's response is simply a shifted impulse response $h(t-\tau)$.



(2) Linearity implies that the response to a scaled impulse input $x(\tau) \delta(t-\tau)$ is $x(\tau) h(t-\tau)$—that is, the impulse response is scaled by the amplitude $x(\tau)$.



Response to Arbitrary Input $x(t)$

The impulse response $h(t)$ characterizes how the system reacts to the simplest possible input: the Dirac delta function $\delta(t)$, which represents a "spike" at a single point in time.


A general input signal $x(t)$ can be decomposed into an infinite sum (or integral) of weighted, shifted impulses:


$$ \begin{align*} x(t) &= \int_{-\infty}^{\infty} x(\tau) \, \delta(t-\tau) \, d\tau \end{align*} $$


Each $\delta(t-\tau)$ in this decomposition contributes a response $h(t-\tau)$.


By superposition (linearity), the system output $y(t)$ is the sum (integral) of these scaled and shifted responses:


$$y(t) = \int_{-\infty}^{\infty} x(\tau) \, h(t-\tau) \, d\tau \; = \; x(t) * h(t)$$



Thus, convolution with $h(t)$ completely characterizes the system’s behavior for any input $x(t)$.


Example:

Let's consider:


$$ \dot{y} + 3 y = 3 x(t) $$


with impulse response:


$$ h(t) = 3e^{-3t}, \quad t \ge 0 $$


If $x(t)$ is a general input (e.g., square wave, ramp, sine wave), we can compute $y(t)$ by:


$$ y(t) = h(t)*x(t) = \int_{0}^{t} 3e^{-3(t-\tau)} x(\tau) d\tau $$


This integral sums the weighted past values of $x(\tau)$, with the weighting determined by how much time has passed (via $3e^{-3(t-\tau)}$).


In [ ]:
% generate a general input

T = 4;

fs = 100;
Ts = 1/fs;
t = 0:Ts:2*T-Ts;

w0 = 2*pi/T;
x = 1/2*square(w0*t);

plot(t,x,'k')
grid on
ylim([-0.6,0.6])
yticks([-0.5,0,0.5])

No description has been provided for this image
In [ ]:
% convolution (circular, will be discussed later)

h = 3*exp(-3*t);
y = cconv(x,h,length(t))*Ts;

plot(t,x,'k'), hold on
plot(t,y,'b', 'linewidth',2), hold off
grid on
ylim([-0.6,0.6])
yticks([-0.5,0,0.5])

No description has been provided for this image




4. Frequency Response to General Input (in Frequency)

4.1. Response to a Sinusoidal Input

A remarkable property of LTI (Linear Time-Invariant) systems is that when the input is a complex exponential of the form


$$ x(t) = e^{j\omega t} $$


the output is the same exponential multiplied by a complex number that depends only on $\omega$:




Let the system's impulse response be $h(t)$. The output is given by convolution:


$$ y(t) = h(t) * x(t) = \int_{-\infty}^{\infty} h(\tau) \, x(t - \tau) \, d\tau $$


Substituting $x(t) = e^{j\omega t}$:


$$ \begin{align*} y(t) &= \int_{-\infty}^{\infty} h(\tau) \, e^{j\omega(t - \tau)} \, d\tau \\\\ &= e^{j\omega t} \underbrace{\int^\infty_{-\infty}h(\tau)e^{-j\omega \tau}d\tau}_{\text{complex function of } \omega} \\\\ &= e^{j\omega t} \cdot H(j\omega) \end{align*} $$


where


$$ H(j\omega) = \int_{-\infty}^{\infty} h(t) e^{-j\omega t} dt $$


is the Fourier Transform of the impulse response $h(t)$.



Definition: Fourier Transform


$$H(j\omega) = \int^\infty_{-\infty}h(\tau)e^{-j\omega \tau}d\tau$$




Interpretation

  • The system does not change the frequency of the input.
  • The output is a scaled and phase-shifted version of the input.

That is,


$$ y(t) = \lvert H(j\omega) \rvert \cdot e^{j\omega t + j \angle H(j\omega)} $$


Here:

  • $\lvert H(j\omega) \rvert$ is the gain (amplitude response)
  • $\angle H(j\omega)$ is the phase shift introduced by the system


$y(t) = H(j\omega)e^{j\omega t}$ rotates with the same angular velocity $\omega$


$$ \begin{align*} y(t) &= H(j\omega) e^{j\omega t}\\ (j\omega) &= \lvert H(j\omega) \rvert e^{j\angle H(j\omega)} \end{align*} $$



Example: Frequency Response of


$$ \dot{y} + 5y = x(t)$$


This is a linear time-invariant (LTI) system. We want to understand how it responds to sinusoidal input using the frequency domain.


We express $H(j\omega)$ in polar form:


$$ \lvert H(j\omega) \rvert = \frac{1}{\sqrt{25 + \omega^2}}, \qquad \angle H(j\omega) = -\tan^{-1}\left( \frac{\omega}{5} \right) $$


In [ ]:
% use lsim

A = -5;
B = 1;
C = 1;
D = 0;
G = ss(A,B,C,D);

w = pi;
t = linspace(0,2*pi,200);

x0 = 0;

x = sin(w*t);

[y,tout] = lsim(G,x,t,x0);

plot(t,x), hold on
plot(tout,y), hold off
grid on
xlabel('t')
axis tight, ylim([-1,1])
leg = legend('input','output');
set(leg, 'fontsize', 12)

No description has been provided for this image
In [ ]:
% sinusoidal inputs with different w

A = -5;
B = 1;
C = 1;
D = 0;
G = ss(A,B,C,D);

x0 = 0;

t = linspace(0,2*pi,200);

W = [1,5,10,20];
for w = W
    x = sin(w*t);
    [y,tout] = lsim(G,x,t,x0);
    plot(tout,y), hold on
end

hold off
grid on
axis tight, ylim([-0.3,0.3])
xlabel('t')
leg = legend('\omega = 1','\omega = 5','\omega = 10','\omega = 20');
set(leg, 'fontsize', 12)

No description has been provided for this image

4.2. Response to a Periodic Input (in the Frequency Domain)


4.2.1. Periodic Signal and Harmonics

A signal is periodic if it satisfies:


$$x(t) = x(t + T) $$


where:


$$ \begin{align*} {\omega_0} &= 2\pi f \quad \text{is the fundamental frequency }\\ T &= \frac{2\pi}{\omega_0} \quad \;\text{is the fundamental period } \end{align*} $$



Periodic signals can be represented as an infinite sum of complex exponentials using the Fourier series.



Fourier series


$$ \begin{align*} x(t) &= a_0 + a_1e^{j\omega_0t} + a_2e^{j2\omega_0t} + a_3e^{j3\omega_0t} + \cdots \\\\ &= \sum\limits^\infty_{k=-\infty}a_ke^{jk\omega_0t} \end{align*} $$


This representation decomposes the signal into harmonic components - sinusoids with frequencies that are integer multiples of the fundamental frequency.



Let $\omega_0 = \frac{2\pi}{T}$ be the fundamental frequency. The complex exponentials $e^{jk\omega_0 t}$ are orthogonal over a period $T$:


$$ \int_t^{t+T} e^{jk_1\omega_0 \tau} \cdot e^{-jk_2\omega_0 \tau} \, d\tau = \begin{cases} 0 & \text{if } k_1 \ne k_2 \\ T & \text{if } k_1 = k_2 \end{cases} \quad \Rightarrow \quad T\delta[k_1 - k_2] $$


This orthogonality allows us to extract the $k$-th coefficient $a_k$ in a Fourier series.


Fourier Series Expansion

Assume $x(t)$ is periodic with period $T$:


$$ x(t) = \sum_{k=-\infty}^{\infty} a_k e^{j k \omega_0 t}, \qquad \omega_0 = \frac{2\pi}{T} $$


Finding the Coefficients $a_k$

Multiply both sides by $e^{-j n \omega_0 t}$ and integrate over one period:


$$ \begin{align*} \int_T x(t)e^{-jn\omega_0t}dt &= \int_T \sum\limits^\infty_{k=-\infty}a_ke^{jk\omega_0t} \cdot e^{-jn\omega_0t}dt \\\\ &= \sum\limits^\infty_{k=-\infty}a_k\int_T e^{j(k-n)\omega_0t}dt \\\\ &= a_nT \end{align*} $$


Thus:


$$ a_n = \frac{1}{T} \int_T x(t)e^{-j n \omega_0 t} dt $$


Summary

  • Analysis (extracting $a_k$):

$$ a_k = \frac{1}{T} \int_T x(t) e^{-j k \omega_0 t} \, dt $$


  • Synthesis (rebuilding $x(t)$):

$$ x(t) = \sum_{k=-\infty}^{\infty} a_k e^{j k \omega_0 t} $$


This framework applies to any periodic signal that meets Dirichlet conditions.



Example: Fourier Series of a Triangle Wave

The triangle wave is a periodic, piecewise linear signal that is odd-symmetric and continuous, often defined over the interval $[-T/2, T/2]$. Due to its symmetry and shape, its Fourier series contains only odd harmonics (i.e., $k = \pm 1, \pm 3, \pm 5, \dots$), and the coefficients decay with $1/k^2$.


The Fourier series expansion in complex exponential form is:


$$ x(t) = \sum_{\substack{k = -\infty \\ k \; \text{odd}}}^{\infty} \frac{-1}{2k^2\pi^2} e^{j 2\pi k t} $$


Alternatively, using the fundamental frequency $\omega_0 = 2\pi$, this becomes:


$$ x(t) = \sum_{\substack{k = -\infty \\ k \; \text{odd}}}^{\infty} \frac{-2}{\omega_0^2 k^2} e^{j \omega_0 k t} $$


Key Features:

  • Odd harmonics only: The triangle wave contains no even harmonics due to its odd symmetry.
  • Spectral decay: Coefficients decay as $1/k^2$, which means the triangle wave has a smooth appearance (unlike square waves, whose coefficients decay as $1/k$).
  • Real-valued signal: Since $x(t)$ is real, the coefficients satisfy $a_{-k} = a_k^*$, ensuring complex exponentials combine into real sinusoids.

In [ ]:
T = 1;

fs = 100;
t = 0:1/fs:2*T-1/fs;

w0 = 2*pi/T;
x = 1/8*sawtooth(w0*t, 1/2);

% xhat = Fourier series of triangle ware
xhat = zeros(size(t));

for k = -3:2:3
    xhat = xhat - 2/(w0^2*k^2)*exp(1j*w0*k*t);
end

plot(t, x, 'k'), hold on
plot(t, xhat, 'b', 'linewidth', 2), hold off
grid on
ylim([-2/8,2/8])
xticks([0,1,2])
yticks([-1/8,0,1/8])
Warning: Imaginary parts of complex X and/or Y arguments ignored

No description has been provided for this image

Example: Fourier Series of a Square Wave

A square wave is a periodic, piecewise constant signal that alternates between two values (e.g., $+1$ and $-1$), often with odd symmetry (i.e., $x(-t) = -x(t)$). It is defined over a fundamental period $T$, with fundamental frequency $\omega_0 = \frac{2\pi}{T}$.


Its Fourier series expansion in complex exponential form is:


$$ x(t) = \sum_{\substack{k = -\infty \\ k \; \text{odd}}}^{\infty} \frac{1}{jk\pi} e^{j 2\pi k t} $$


or equivalently using $\omega_0 = 2\pi$:


$$ x(t) = \sum_{\substack{k = -\infty \\ k \; \text{odd}}}^{\infty} \frac{2}{j\omega_0 k} e^{j \omega_0 k t} $$


Key Characteristics:

  • Odd harmonics only: Square waves only contain odd multiples of the fundamental frequency due to odd symmetry and half-wave symmetry.
  • Spectral decay: Coefficients decay as $1/k$, which means the signal has sharp transitions (discontinuities).
  • Imaginary denominator: Indicates that when converted to sine-cosine form, the square wave becomes a sum of sines, since the imaginary parts of complex exponentials give sine functions.
  • Real signal: Despite using complex exponentials, the result is real-valued due to the symmetry $a_{-k} = a_k^*$.

In [ ]:
T = 1;

fs = 100;
t = 0:1/fs:2*T-1/fs;

w0 = 2*pi/T;
x = 1/2*square(w0*t);

% xhat = Fourier series of square ware
xhat = zeros(size(t));

for k = -5:2:5
    xhat = xhat + 2/(1j*w0*k)*exp(1j*w0*k*t);
end

plot(t, x, 'k'), hold on
plot(t, xhat, 'b', 'linewidth', 2), hold off
grid on
ylim([-5/8,5/8])
xticks([0,1,2])
yticks([-1/2,0,1/2])
Warning: Imaginary parts of complex X and/or Y arguments ignored

No description has been provided for this image

4.2.2. LTI System Response to Harmonic Inputs



One of the most powerful properties of LTI systems is that complex exponentials are eigenfunctions. That is, for an input:


$$ x_k(t) = e^{j\omega_k t} $$


the output is:


$$ y_k(t) = H(j\omega_k) e^{j\omega_k t} $$


where $H(j\omega_k)$ is the system’s frequency response evaluated at frequency $\omega_k$.


Thus, for a general periodic input:


$$ x(t) = \sum_{k=-\infty}^{\infty} a_k e^{j\omega_0 k t} $$


the output is simply:


$$ y(t) = \sum_{k=-\infty}^{\infty} a_k H(j\omega_0 k) e^{j\omega_0 k t} $$


This is a filtered version of the input: each harmonic is scaled and phase-shifted according to the system's frequency response at that harmonic.



Worked Example: First-Order System

Consider the first-order LTI system:


$$ \dot{y} + \frac{1}{\tau} y = \frac{1}{\tau} x(t) $$


Let the input be a complex exponential:


$$ x(t) = e^{j\omega t} $$


We look for a steady-state solution of the form:


$$ y(t) = A e^{j(\omega t + \phi)} $$


Substituting into the ODE:


$$ \frac{d}{dt} \left(A e^{j(\omega t + \phi)}\right) + \frac{1}{\tau} A e^{j(\omega t + \phi)} = \frac{1}{\tau} e^{j\omega t} $$


yields:


$$ \left(j\omega + \frac{1}{\tau} \right) A e^{j\phi} = \frac{1}{\tau} $$


Solving for the amplitude and phase:


$$ A = |H(j\omega)| = \frac{1}{\sqrt{(\tau \omega)^2 + 1}}, \quad \phi = -\tan^{-1}(\tau \omega) $$


Thus, the frequency response is:


$$ H(j\omega) = \frac{1}{1 + j\tau\omega} $$


with magnitude and phase:


$$ |H(j\omega)| = \frac{1}{\sqrt{1 + (\tau\omega)^2}}, \qquad \angle H(j\omega) = -\tan^{-1}(\tau \omega) $$


Interpretation

  • Each frequency component of a periodic signal is attenuated and phase-shifted by the system's response at that frequency.
  • The output spectrum is a modulated version of the input spectrum.
  • The system behaves like a filter, shaping the harmonic content of the signal.









In [ ]:
T = 1;

fs = 100;
t = 0:1/fs:2*T-1/fs;

w0 = 2*pi/T;
x = 1/2*square(w0*t);

% xhat = Fourier series of square ware
xhat = zeros(size(t));

for k = -19:2:19
    xhat = xhat + 2/(1j*w0*k)*exp(1j*w0*k*t);
end

plot(t, x, 'k'), hold on
plot(t, xhat, 'b', 'linewidth', 2), hold off
grid on
ylim([-5/8,5/8])
xticks([0,1,2])
yticks([-1/2,0,1/2])
Warning: Imaginary parts of complex X and/or Y arguments ignored

No description has been provided for this image
In [ ]:
tau = 1/5;

yhat = zeros(size(t));

for k = -19:2:19
    w = w0*k;
    A = 1./abs(1j*tau*w + 1);
    Ph = -phase(1j*tau*w + 1);
    yhat = yhat + A*2/(1j*w0*k)*exp(1j*(w*t + Ph));
end

plot(t, x, 'k'), hold on
plot(t, yhat, 'g', 'linewidth', 2), hold off
grid on
ylim([-5/8,5/8])
xticks([0,1,2])
yticks([-1/2,0,1/2])
leg = legend('square', 'output');
set(leg, 'fontsize', 12, 'location', 'best')
Warning: Imaginary parts of complex X and/or Y arguments ignored

No description has been provided for this image

4.3. Response to a General Input (Aperiodic Signal) in Frequency Domain

From Periodic to Aperiodic

An aperiodic signal can be interpreted as a periodic signal with infinite period. That is:


$$x_T(t) = \sum\limits^\infty_{k=-\infty}x(t+kT)$$


$$x(t) = \lim\limits_{T\rightarrow \infty}x_T(t)$$



This periodic extension allows the use of Fourier series, which then transition to the Fourier transform as $T \to \infty$.


Example: Periodic Square Wave

Let a aperiodic signal be defined as:


$$ x(t) = \begin{cases} 1 & |t| < S \\ 0 & S < |t| < \frac{T}{2} \end{cases} $$



Then periodic extension to periodic square wave $x_T(t)$


$$x_T(t) = \sum\limits^\infty_{k=-\infty}x(t+kT)$$



Then compute the Fourier series coefficients:


$$ \begin{align*} \omega_0 &= {2\pi \over T}\\\\ a_k &= {1 \over T} \int^{T \over 2}_{-T \over 2}x_T(t)e^{-j\omega_0kt}dt \\ &= {1 \over T} \bigg[\int^s_{-s}1 \cdot e^{-j\omega_0kt}dt \bigg] = {1 \over T}\bigg[-{e^{-j\omega_0kt} \over j\omega_0k} \bigg\rvert^s_{-s} \bigg] \\ &= -{1 \over T} \bigg[ {e^{-j\omega_0ks} \over j\omega_0k} - {e^{j\omega_0ks} \over j\omega_0k} \bigg] = {1 \over T} \bigg[{-2 \over j\omega_0k} \cdot {{e^{-j\omega_0ks} - e^{j\omega_0ks}} \over 2} \bigg] \\ &= {1 \over T} {-2 \over j \omega_0 k}(-j \sin (\omega_0 k s)) = {1 \over T}{2\sin (\omega_0 k s) \over \omega_0 k} \\ &= {2 \over T} {\sin(\omega_0 k s) \over \omega_0k} \end{align*} $$


Multiply by $T$ to isolate the spectral content:


$$ T a_k = \frac{2 \sin(\omega_0 k s)}{\omega_0 k} $$


Let $\omega = \omega_0 k$, then:


$$ T a_k = \frac{2 \sin(\omega s)}{\omega} $$



Limit as $T \to \infty$

As the period $T$ becomes infinite:

  • The fundamental frequency $\omega_0 = \frac{2\pi}{T} \to 0$
  • The discrete set of harmonic frequencies becomes dense, forming a continuum
  • The scaled coefficients $T a_k$ approach the continuous Fourier transform:

$$ \lim_{T \to \infty} T a_k = X(j\omega) = \int_{-\infty}^\infty x(t) e^{-j\omega t} dt $$


For the square pulse:


$$ X(j\omega) = \frac{2 \sin(\omega s)}{\omega} $$


This is the well-known sinc-shaped spectrum of a rectangular pulse.


Summary

  • Periodic signal → sum of discrete harmonics via Fourier series
  • Aperiodic signal → continuous frequency spectrum via Fourier transform
  • As $T \to \infty$, the discrete coefficients $a_k$ become samples of an envelope function:

$$ T a_k \to X(j\omega) = \frac{2 \sin(\omega s)}{\omega} $$


$$Ta_k = {2\text{sin} \omega s \over \omega} \bigg\rvert_{\omega=k\omega_0}$$


  • That is, with $\omega$ thought of as a continuous variable, the set of Fourier series coefficients approaches the envelop function as $T \rightarrow \infty$

  • Aperiodic signal has all the frequency components instead of discrete harmonic components

  • The Fourier series becomes a Fourier transform:


$$ x(t) = \sum_k a_k e^{j \omega_0 k t} \quad \to \quad x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j \omega t} d\omega $$


This forms the foundation for analyzing arbitrary signals and their behavior in LTI systems using frequency response.





5. Fourier Transform

From Periodic to Aperiodic

As the period $T \to \infty$, the discrete frequencies $\omega_0 = \frac{2\pi}{T}$ become infinitely close, turning the summation into an integral:


$$ \begin{align*} x(t) &= \sum_{k=-\infty}^{\infty} a_k e^{j\omega_0 k t} \\\\ &= \sum_{k=-\infty}^{\infty} \frac{X(j\omega_0 k)}{T} e^{j\omega_0 k t} \\\\ &= \frac{1}{2\pi} \sum_{k=-\infty}^{\infty} X(j\omega_0 k) e^{j\omega_0 k t} \, \omega_0 \\\\ &\to \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} d\omega \quad \text{as } \; \;T \to \infty \end{align*} $$



Fourier Transform Pair

This leads to the Fourier Transform, which maps time-domain signals into the frequency domain:


Analysis (forward transform):


$$ X(j\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} dt $$


Synthesis (inverse transform):


$$ x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} d\omega $$


Response of LTI Systems in Frequency Domain

Given an LTI system with impulse response $h(t)$ and input $x(t)$:

(1) Time-domain:


$$y(t) = x(t) * h(t)$$


(2) Frequency-domain:


$$ \mathscr{F}\{y(t)\} = \mathscr{F}\{x(t)*h(t)\} = X(j\omega)H(j\omega) $$


So,


$$ y(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega)H(j\omega) e^{j\omega t} d\omega $$


Summary

  • Fourier Series: for periodic signals
  • Fourier Transform: for aperiodic signals
  • LTI system response is multiplication in frequency domain:

$$ y(t) = x(t) * h(t) \quad \Longleftrightarrow \quad Y(j\omega) = X(j\omega) H(j\omega) $$


This is a cornerstone of signal processing and control theory.


5.1. Impulse Response

The impulse response $h(t)$ fully characterizes a linear time-invariant (LTI) system. It is the output of the system when the input is a Dirac delta function $\delta(t)$:


$$ \delta(t) \quad \mathop\longrightarrow^{\text{LTI}} \quad h(t) $$


Fourier Transform of the Delta Function

The Dirac delta function has a very special Fourier transform:


$$ \int_{-\infty}^{\infty} \delta(t) e^{-j\omega t} dt = 1 $$


So,


$$ \delta(t) \quad \longleftrightarrow \quad 1 \quad \text{(in frequency domain)} $$



This means $\delta(t)$ contains all frequencies equally, with unit magnitude and zero phase shift.


Why Impulse Response is Powerful

  • Excitation by all frequencies: Since $\delta(t)$ contains all frequencies, applying it to an LTI system "probes" the system's response at every frequency.

  • Impulse basically excites a system with all the frequency of $𝑒^{𝑗\omega 𝑡}$

  • Time domain: $h(t)$ describes how the system reacts over time.

  • Frequency domain: Its Fourier transform $H(j\omega)$ tells how each frequency is:

    • Scaled (magnitude response $\lvert H(j\omega) \rvert$)
    • Shifted (phase response $\angle H(j\omega)$)


Summary

Impulse response is the most fundamental descriptor of an LTI system both in time and frequency.

  • Impulse response $h(t)$ $\leftrightarrow$ frequency response $H(j\omega)$
  • $H(j\omega)$ tells us how the system filters sinusoidal components of any input
  • Impulse response contains the information on how much magnitude and phase are filtered via the LTI system at all the frequency
  • The convolution $y(t) = h(t) * x(t)$ in time is equivalent to multiplication in frequency:

$$ x(t) * h(t) \quad \longleftrightarrow \quad X(j\omega) \cdot H(j\omega) $$




5.2. Frequency Response (Frequency Sweep)

A frequency sweep is another way to characterize an LTI system - equivalent in information to the impulse response, but presented in the frequency domain.


Suppose we input a complex exponential at a single frequency:


$$ x(t) = e^{j\omega t} $$


and consider the LTI system:


$$ \dot{y} + 5y = 5x(t) $$


Substituting the input:


$$ \dot{y} + 5y = 5e^{j\omega t} $$


Assume the output has the same frequency as the input (because LTI systems do not create new frequencies):


$$ y(t) = A e^{j(\omega t + \phi)} = A e^{j\omega t} e^{j\phi} $$


Take the derivative:


$$ \dot{y}(t) = j\omega A e^{j\omega t} e^{j\phi} $$


Substitute into the differential equation:


$$ (j\omega + 5) A e^{j\omega t} e^{j\phi} = 5 e^{j\omega t} $$


Cancel $e^{j\omega t}$ from both sides:


$$ (j\omega + 5) A e^{j\phi} = 5 $$


Take the magnitude:


$$ A = \frac{5}{|j\omega + 5|} = \frac{5}{\sqrt{\omega^2 + 25}} $$


Take the phase:


$$ \phi = -\angle(j\omega + 5) = -\tan^{-1}\left( \frac{\omega}{5} \right) $$


Frequency Response

The frequency response of the system is:


$$ H(j\omega) = \frac{5}{j\omega + 5} $$


  • Magnitude response:

$$ |H(j\omega)| = \frac{5}{\sqrt{\omega^2 + 25}} $$


  • Phase response:

$$ \angle H(j\omega) = -\tan^{-1} \left( \frac{\omega}{5} \right) $$


Summary

By sweeping over $\omega$ and observing the amplitude and phase of the output, we obtain:

  • The gain at each frequency: how much the system amplifies or attenuates each frequency
  • The phase shift at each frequency: how the output is delayed or advanced relative to the input

This is precisely the frequency response $H(j\omega)$ of the LTI system.


In [ ]:
w = 0.1:0.1:100;

A = 5./abs(1j*w+5);
P = -angle(1j*w+5)*180/pi;

subplot(2,1,1), plot(w, A, 'linewidth', 2)
yticks([0,0.5,1])
ylabel('|H(j\omega)|', 'fontsize', 13)
xlabel('\omega', 'fontsize', 15)

subplot(2,1,2), plot(w, P, 'linewidth', 2)
yticks([-90, -45, 0])
ylabel('\angle H(j\omega)', 'fontsize', 13)
xlabel('\omega', 'fontsize', 15)

% later, we will see that this is kind of a bode plot

No description has been provided for this image



5.3. Frequency Response of a Second-Order System

Consider the second-order ODE:


$$ \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2 y = \omega_n^2 \; x(t) $$


We input a complex exponential:


$$ x(t) = e^{j\Omega t} $$


The differential equation becomes:


$$ \ddot{y} + 2\zeta\omega_n\dot{y} + \omega_n^2 y = \omega_n^2 \; e^{j\Omega t} $$


(1) Step 1: Assume the Steady-State Output

Because $x(t)$ is a complex exponential and LTI systems preserve frequency:


$$ y(t) = A e^{j(\Omega t + \phi)} = A e^{j\phi} e^{j\Omega t} $$


(2) Step 2: Substitute into the Differential Equation

Take derivatives:


$$ \dot{y} = j\Omega A e^{j\phi} e^{j\Omega t}, \qquad \ddot{y} = -\Omega^2 A e^{j\phi} e^{j\Omega t} $$


Substitute into the ODE:


$$ \left( -\Omega^2 + j 2\zeta \omega_n \Omega + \omega_n^2 \right) A e^{j\phi} e^{j\Omega t} = \omega_n^2 e^{j\Omega t} $$


Cancel $e^{j\Omega t}$:


$$ A e^{j\phi} = \frac{\omega_n^2}{-\Omega^2 + j 2\zeta \omega_n \Omega + \omega_n^2} = \frac{1}{1-\left(\frac{\Omega}{\omega_n}\right)^2 + j 2\zeta \left(\frac{\Omega}{\omega_n}\right)} $$


(3) Step 3: Express in Terms of Frequency Ratio $\gamma = \frac{\Omega}{\omega_n}$

The magnitude (gain) becomes:


$$ A = \left| \frac{1}{1 - \left( \frac{\Omega}{\omega_n} \right)^2 + j 2\zeta \left( \frac{\Omega}{\omega_n} \right)} \right| = \frac{1}{ \sqrt{ \left(1 - \gamma^2 \right)^2 + 4\zeta^2 \gamma^2 } } $$


The phase shift is:


$$ \phi = -\tan^{-1} \left( \frac{2\zeta \frac{\Omega}{\omega_n}}{1-\left(\frac{\Omega}{\omega_n}\right)^2} \right) =-\tan^{-1} \left( \frac{2\zeta \gamma}{1-\gamma^2} \right) $$


Summary

For the second-order system with input $x(t) = e^{j\Omega t}$:

  • Frequency ratio: $\gamma = \frac{\Omega}{\omega_n}$
  • Amplitude response (gain):

$$ |H(j\Omega)| = \frac{1}{ \sqrt{ (1 - \gamma^2)^2 + 4\zeta^2 \gamma^2 } } $$


  • Phase response:

$$ \angle H(j\Omega) = -\tan^{-1} \left( \frac{2\zeta \gamma}{1 - \gamma^2} \right) $$


This is the standard frequency response of a damped second-order system.


In [ ]:
r = 0:0.05:4;

zeta = 0.1:0.2:1;
A = [];
for i = 1:length(zeta)
    A(i,:) = 1./sqrt((1-r.^2).^2 + (2*zeta(i)*r).^2);
end

plot(r, A)
xlabel('\gamma')
ylabel('A')
legend('0.1','0.3','0.5','0.7','0.9')

No description has been provided for this image
In [ ]:
phi = [];
for i = 1:length(zeta)
    phi(i,:) = -atan2((2*zeta(i).*r),(1-r.^2));
end

plot(r,phi*180/pi)
xlabel('\gamma')
ylabel('\phi')
legend('0.1','0.3','0.5','0.7','0.9')

No description has been provided for this image



5.4. Resonance in Second-Order Systems

  • When the input frequency $\Omega$ approaches a specific value, the amplitude of the output response becomes maximal.
  • This phenomenon is called resonance.

What is Resonance?

The resonance frequency $\Omega_r$ is the input frequency that maximizes the amplitude of the steady-state response:


$$ A = \left| H(j\Omega) \right| = \frac{1}{ \sqrt{ (1 - \gamma^2)^2 + 4\zeta^2 \gamma^2 } }, \quad \gamma = \frac{\Omega}{\omega_n} $$


To find the resonance frequency, minimize the denominator:


$$ (1 - \gamma^2)^2 + 4\zeta^2 \gamma^2 $$


Taking derivative with respect to $\gamma$ and solving yields:


$$ \gamma_r = \sqrt{1 - 2\zeta^2} \quad \Longrightarrow \quad \Omega_r = \omega_n \sqrt{1 - 2\zeta^2} $$


Key Points

  • Resonance frequency:

$$ \Omega_r = \omega_n \sqrt{1 - 2\zeta^2} $$


  • Natural frequency:

$$ \omega_n $$


  • For small damping ($\zeta \ll 1$), $\Omega_r \approx \omega_n$

  • Resonance frequency is generally different from natural frequency, but they often are close enough


Physical Meaning

  • The system amplifies input signals near $\Omega_r$.
  • If damping $\zeta$ is low, this amplification can be very large, leading to large oscillations.

In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('CzJ8EKi11pU', width = "560", height = "315")
Out[ ]:
In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('Mqt_Amkw-0Y', width = "560", height = "315")
Out[ ]:
In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('TYA8C85CdCc', width = "560", height = "315")
Out[ ]:
In [1]:
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')