Example: Time-Varying Frequency Response Maps with sid.freq_map¶
sid.freq_map estimates how a plant's frequency response evolves
over time by running spectral analysis on overlapping segments. It
is the right tool for systems whose dynamics drift (thermally
softening structures, aging actuators, fuel-mass change in
rockets), snap (latch release, damage events), or depend on
operating point (nonlinearities like Duffing).
This notebook works through four scenarios on physical SMD plants:
- LTI baseline — constant 2-mass chain; the map should look time-invariant.
- Continuous LTV — the first spring's stiffness ramps linearly in time, simulating thermal softening.
- Discrete LTV — the first spring's stiffness snaps at
t = T/2, simulating a sudden structural change. - Duffing hardening — a nonlinear SDOF whose effective resonance frequency rises with amplitude. A ramped-amplitude excitation causes the apparent frequency to drift upward over the record.
Translated from exampleFreqMap.m.
import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd, util_msd_ltv, util_msd_nl
import sid
%matplotlib inline
1. LTI baseline: constant 2-mass chain¶
Plant m = [1,1], k = [100, 80], c = [2, 2]. Force is applied
at mass 1 and we measure mass 1's position. The frequency map
should be constant along the time axis — the test of whether
segment-based analysis produces a stable picture in the absence
of any drift.
rng = np.random.default_rng(10)
m = np.array([1.0, 1.0])
k = np.array([100.0, 80.0])
c = np.array([2.0, 2.0])
F = np.array([[1.0], [0.0]])
Ts = 0.01
N = 4000
Ad, Bd = util_msd(m, k, c, F, Ts)
u = rng.standard_normal(N)
x = np.zeros((N + 1, 4))
for step in range(N):
x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step]
y = x[1:, 0] + 5e-4 * rng.standard_normal(N)
result_lti = sid.freq_map(y, u, segment_length=512, overlap=384, sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
sid.map_plot(result_lti, plot_type='magnitude', ax=ax)
ax.set_title('LTI 2-mass chain: magnitude should be constant along time')
plt.tight_layout()
plt.show()
2. Continuous LTV: ramping first-spring stiffness¶
Now k₁(t) ramps linearly from 200 N/m down to 20 N/m over
the record — a 10× softening. The first mode drops from roughly
7.2 rad/s to 3.1 rad/s, the second from 17.6 rad/s to
13.1 rad/s. We build the per-step (Ad, Bd) stack with
util_msd_ltv (pass the stiffness as a (2, N) array) and
simulate the LTV state recursion.
rng2 = np.random.default_rng(20)
k_tv = np.zeros((2, N))
k_tv[0, :] = np.linspace(200.0, 20.0, N) # soften mass 1's spring
k_tv[1, :] = 80.0 # k₂ stays constant
m_tv = np.broadcast_to(m[:, None], (2, N))
c_tv = np.broadcast_to(c[:, None], (2, N))
Ad_tv, Bd_tv = util_msd_ltv(m_tv, k_tv, c_tv, F, Ts) # (4,4,N), (4,1,N)
u_tv = rng2.standard_normal(N)
x_tv = np.zeros((N + 1, 4))
for step in range(N):
x_tv[step + 1] = Ad_tv[:, :, step] @ x_tv[step] + Bd_tv[:, 0, step] * u_tv[step]
y_tv = x_tv[1:, 0] + 5e-4 * rng2.standard_normal(N)
result_ramp = sid.freq_map(y_tv, u_tv, segment_length=256, overlap=192, sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
sid.map_plot(result_ramp, plot_type='magnitude', ax=ax)
ax.set_title('Continuous LTV: k₁(t) ramps 200 → 20 (first mode drifts 7 → 3 rad/s)')
plt.tight_layout()
plt.show()
3. Discrete LTV: step change in stiffness¶
Same plant but k₁ snaps from 200 N/m to 40 N/m at the middle
of the record. This is a caricature of a sudden structural event
(cable cut, latch release, damage). The frequency map should show
an abrupt transition at t = T/2; the smoothness of the
transition depends on the segment length — longer segments straddle
the jump and smear it across several time bins.
rng3 = np.random.default_rng(30)
k_step = np.zeros((2, N))
k_step[0, :N // 2] = 200.0 # stiff first half
k_step[0, N // 2:] = 40.0 # soft second half
k_step[1, :] = 80.0
Ad_step, Bd_step = util_msd_ltv(m_tv, k_step, c_tv, F, Ts)
u_step = rng3.standard_normal(N)
x_step = np.zeros((N + 1, 4))
for step in range(N):
x_step[step + 1] = Ad_step[:, :, step] @ x_step[step] + Bd_step[:, 0, step] * u_step[step]
y_step = x_step[1:, 0] + 5e-4 * rng3.standard_normal(N)
result_step = sid.freq_map(y_step, u_step, segment_length=256, overlap=192, sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
sid.map_plot(result_step, plot_type='magnitude', ax=ax)
ax.set_title('Discrete LTV: step change k₁: 200 → 40 at t = T/2')
plt.tight_layout()
plt.show()
Coherence map¶
Coherence shows how the signal-to-noise ratio evolves over time. For the ramping-stiffness run, coherence tends to drop as the plant softens and the resonance slides toward DC where less excitation reaches.
fig, ax = plt.subplots(figsize=(8, 5))
sid.map_plot(result_ramp, plot_type='coherence', ax=ax)
ax.set_title('Coherence map: ramping stiffness')
plt.tight_layout()
plt.show()
BT vs Welch algorithm¶
sid.freq_map supports two algorithms:
'bt'(Blackman-Tukey, default): correlogram-based, configurable lag window.'welch': Welch averaged periodogram, sub-segment FFT averaging.
Both produce similar magnitudes on this plant; BT is smoother (fewer bins) while Welch gives finer frequency resolution at the cost of higher variance per bin.
result_bt = sid.freq_map(y_tv, u_tv, segment_length=256, algorithm='bt', sample_time=Ts)
result_welch = sid.freq_map(y_tv, u_tv, segment_length=256, algorithm='welch', sample_time=Ts)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sid.map_plot(result_bt, plot_type='magnitude', ax=axes[0])
axes[0].set_title('Blackman-Tukey')
sid.map_plot(result_welch, plot_type='magnitude', ax=axes[1])
axes[1].set_title('Welch')
plt.tight_layout()
plt.show()
Segment length and overlap tuning¶
Shorter segments give better time resolution but worse frequency resolution. More overlap gives a smoother time axis but requires more computation. On the ramping LTV plant, the short-segment view captures the smooth drift more clearly, while the long-segment view averages several time windows and shows a blurrier transition.
result_short = sid.freq_map(y_tv, u_tv, segment_length=128, overlap=96, sample_time=Ts)
result_long = sid.freq_map(y_tv, u_tv, segment_length=512, overlap=384, sample_time=Ts)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
sid.map_plot(result_short, plot_type='magnitude', ax=axes[0])
axes[0].set_title('Short segments (L = 128): good time resolution')
sid.map_plot(result_long, plot_type='magnitude', ax=axes[1])
axes[1].set_title('Long segments (L = 512): good frequency resolution')
plt.tight_layout()
plt.show()
Time-series mode: evolving output spectrum¶
With u = None, sid.freq_map estimates how the output spectrum
changes over time. We re-run the ramping-stiffness plant with
fresh noise and hand only the position record to freq_map.
rng4 = np.random.default_rng(40)
u_nu = rng4.standard_normal(N)
x_nu = np.zeros((N + 1, 4))
for step in range(N):
x_nu[step + 1] = Ad_tv[:, :, step] @ x_nu[step] + Bd_tv[:, 0, step] * u_nu[step]
y_nu = x_nu[1:, 0]
result_ts = sid.freq_map(y_nu, None, segment_length=256, overlap=192, sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
sid.map_plot(result_ts, plot_type='spectrum', ax=ax)
ax.set_title('Time-series: output spectrum drifts as k₁ softens')
plt.tight_layout()
plt.show()
4. Duffing hardening oscillator¶
A nonlinear SDOF with cubic stiffness:
m·ẍ + c·ẋ + k_lin·x + k_cub·x³ = u(t)For a hardening spring (k_cub > 0), the effective resonance
frequency rises with amplitude:
ω_eff(x) ≈ sqrt((k_lin + 3·k_cub·x²) / m)We excite the oscillator with a ramped-amplitude white force (amplitude grows linearly from 0.5 to 10 over the record). As the typical displacement grows, the apparent resonance frequency should drift upward — a signature purely of the nonlinearity, not of any explicit time variation in the plant parameters.
This uses util_msd_nl, the RK4-based nonlinear simulator in
util_msd.py.
rng5 = np.random.default_rng(50)
m_nl = np.array([1.0])
k_lin = np.array([100.0])
k_cubic = np.array([5e4]) # strong hardening
c_nl = np.array([2.0])
F_nl = np.array([[1.0]])
amp = np.linspace(0.5, 10.0, N)
u_nl = (amp * rng5.standard_normal(N)).reshape(-1, 1)
x_nl = util_msd_nl(m_nl, k_lin, k_cubic, c_nl, F_nl, Ts, u_nl, substeps=4)
y_nl = x_nl[1:, 0]
print(f'Typical early amplitude: {np.std(y_nl[:N//4]):.4f} m')
print(f'Typical late amplitude: {np.std(y_nl[-N//4:]):.4f} m')
print(f'Linearized ω_n (small amplitude): '
f'{np.sqrt(k_lin[0] / m_nl[0]):.2f} rad/s')
print(f'Effective ω_n at 0.05 m amplitude: '
f'{np.sqrt((k_lin[0] + 3 * k_cubic[0] * 0.05**2) / m_nl[0]):.2f} rad/s')
result_nl = sid.freq_map(y_nl, u_nl[:, 0], segment_length=256, overlap=192,
sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
sid.map_plot(result_nl, plot_type='magnitude', ax=ax)
ax.set_title('Duffing hardening: apparent resonance rises with amplitude')
plt.tight_layout()
plt.show()
Typical early amplitude: 0.0083 m Typical late amplitude: 0.0315 m Linearized ω_n (small amplitude): 10.00 rad/s Effective ω_n at 0.05 m amplitude: 21.79 rad/s
