Spectrum Analysis

In this example we demonstrate frequency domain analysis using PathSim’s spectrum block. We’ll examine the frequency response of a Butterworth lowpass filter by analyzing its response to a Gaussian pulse.

You can also find this example as a single file in the GitHub repository.

The Spectrum block computes the frequency domain representation of signals during simulation, allowing us to recover the frequency response of linear systems without explicit frequency domain analysis.

First let’s import the Simulation and Connection classes along with the required blocks:

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection
from pathsim.blocks import (
    Scope,
    Spectrum,
    ButterworthLowpassFilter,
    GaussianPulseSource
)
from pathsim.solvers import RKCK54

System Setup

We’ll create a simple signal chain:

1. A Gaussian pulse source with controllable bandwidth

2. A Butterworth lowpass filter

3. Scope and Spectrum blocks to observe time and frequency domain behavior

# Corner frequency of the filter
f = 5

# Blocks that define the system
Src = GaussianPulseSource(f_max=5*f, tau=0.3)
FLT = ButterworthLowpassFilter(f, 10)
Sco = Scope(labels=["pulse", "filtered"])
Spc = Spectrum(freq=np.linspace(0, 2*f, 100),
               labels=["pulse", "filtered"])

blocks = [Src, FLT, Sco, Spc]

The Spectrum block is configured with a frequency array that defines which frequencies to analyze. The GaussianPulseSource provides a rich frequency content suitable for system identification.

# The connections between the blocks
connections = [
    Connection(Src, FLT, Sco, Spc),
    Connection(FLT, Sco[1], Spc[1])
]

We initialize the simulation with the RKCK54 solver (Runge-Kutta-Cash-Karp 5(4) method) for accurate integration of the filter dynamics:

# Initialize simulation
Sim = Simulation(
    blocks,
    connections,
    dt=1e-2,
    log=True,
    Solver=RKCK54,
    tolerance_lte_rel=1e-5,
    tolerance_lte_abs=1e-7
)

07:02:27 - INFO - LOGGING (log: True)
07:02:27 - INFO - BLOCKS (total: 4, dynamic: 2, static: 2, eventful: 0)
07:02:27 - INFO - GRAPH (nodes: 4, edges: 5, alg. depth: 1, loop depth: 0, runtime: 0.079ms)

Now let’s run the simulation:

# Run the simulation for some time
Sim.run(3)

# Plot the time domain results
Sco.plot()

# Plot the frequency domain results
Spc.plot()

plt.show()

07:02:27 - INFO - STARTING -> TRANSIENT (Duration: 3.00s)
07:02:27 - INFO - --------------------   3% | 0.0s<0.1s | 679.0 it/s
07:02:27 - INFO - ####----------------  20% | 0.1s<0.2s | 1444.9 it/s
07:02:27 - INFO - ########------------  40% | 0.1s<0.1s | 1426.9 it/s
07:02:27 - INFO - ############--------  60% | 0.1s<0.0s | 1503.8 it/s
07:02:27 - INFO - ################----  80% | 0.1s<0.0s | 1403.7 it/s
07:02:27 - INFO - #################### 100% | 0.2s<--:-- | 1414.6 it/s
07:02:27 - INFO - FINISHED -> TRANSIENT (total steps: 204, successful: 174, runtime: 153.48 ms)

Frequency Response Recovery

One powerful feature is that we can recover the frequency response of the filter by taking the ratio of the output spectrum to the input spectrum:

# Recover frequency response from spectrum block
freq, (G_pulse, G_filt) = Spc.read()
H_filt_sim = G_filt / G_pulse

We can compare this recovered response with the ideal frequency response calculated from the filter’s state-space representation. The ButterworthLowpassFilter stores its state-space matrices as A, B, C, D:

# Ideal frequency response from filter state-space model
def H(s):
    return np.dot(FLT.C, np.linalg.solve((s*np.eye(FLT.n) - FLT.A), FLT.B)) + FLT.D

H_filt_ideal = np.array([H(2j*np.pi*f) for f in freq]).flatten()

Now let’s plot the comparison between the recovered and ideal frequency responses:

# Plot the frequency domain responses (ideal and recovered)
fig, ax = plt.subplots(nrows=2, tight_layout=True, dpi=120)

ax[0].plot(freq, abs(H_filt_sim), label="recovered")
ax[0].plot(freq, abs(H_filt_ideal), "--", label="ideal")
ax[0].set_xlabel("freq [Hz]")
ax[0].set_ylabel("magnitude")
ax[0].legend()
ax[0].grid(True)

ax[1].plot(freq, np.angle(H_filt_sim), label="recovered")
ax[1].plot(freq, np.angle(H_filt_ideal), "--", label="ideal")
ax[1].set_xlabel("freq [Hz]")
ax[1].set_ylabel("phase [rad]")
ax[1].legend()
ax[1].grid(True)

plt.show()