Example: SISO Frequency Response Estimation¶

This example identifies the frequency response of a physical single-degree-of-freedom spring-mass-damper oscillator from a noisy input/output record using sid.freq_bt (Blackman-Tukey estimator).

Plant. One mass (m = 1 kg) attached to a wall through a spring (k = 100 N/m) and damper (c = 2 N·s/m). The natural frequency is ω_n = sqrt(k/m) = 10 rad/s (≈ 1.59 Hz) and the damping ratio is ζ = c / (2·sqrt(k·m)) = 0.1 (quality factor Q = 5). The input u(t) is a force applied to the mass and the measured output is its position x₁. The plant's static gain from force to position is 1/k = 0.01 m/N, so typical displacements under unit-variance random forcing are a few millimetres.

Translated from exampleSISO.m.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd
import sid

%matplotlib inline

import numpy as np import matplotlib.pyplot as plt from util_msd import util_msd import sid %matplotlib inline

Generate test data¶

We build the continuous-time state-space model of the SMD and discretize it with exact zero-order hold via util_msd. The state vector is [position; velocity]; we measure the position x₁.

In [2]:

rng = np.random.default_rng(42)

# ---- Physical plant: 1-DoF SMD (ω_n = 10 rad/s, ζ = 0.1) ----
m  = np.array([1.0])     # kg
k  = np.array([100.0])   # N/m
c  = np.array([2.0])     # N·s/m
F  = np.array([[1.0]])   # force applied at mass 1
Ts = 0.01                # s (fs = 100 Hz, Nyquist = 50 Hz)
N  = 2048                # number of samples

Ad, Bd = util_msd(m, k, c, F, Ts)   # Ad: (2, 2), Bd: (2, 1)

# ---- Simulate the plant: x[k+1] = Ad x[k] + Bd u[k] ----
u = rng.standard_normal(N)            # white-force excitation
x = np.zeros((N + 1, 2))              # state trajectory [pos; vel]
for step in range(N):
    x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step]

y_clean = x[1:, 0]                    # measured position x₁ (m)
y = y_clean + 2e-4 * rng.standard_normal(N)   # + sensor noise (~0.2 mm)

rng = np.random.default_rng(42) # ---- Physical plant: 1-DoF SMD (ω_n = 10 rad/s, ζ = 0.1) ---- m = np.array([1.0]) # kg k = np.array([100.0]) # N/m c = np.array([2.0]) # N·s/m F = np.array([[1.0]]) # force applied at mass 1 Ts = 0.01 # s (fs = 100 Hz, Nyquist = 50 Hz) N = 2048 # number of samples Ad, Bd = util_msd(m, k, c, F, Ts) # Ad: (2, 2), Bd: (2, 1) # ---- Simulate the plant: x[k+1] = Ad x[k] + Bd u[k] ---- u = rng.standard_normal(N) # white-force excitation x = np.zeros((N + 1, 2)) # state trajectory [pos; vel] for step in range(N): x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step] y_clean = x[1:, 0] # measured position x₁ (m) y = y_clean + 2e-4 * rng.standard_normal(N) # + sensor noise (~0.2 mm)

Estimate frequency response using Blackman-Tukey¶

For a lightly damped structure (Q = 5), the resonance bandwidth is narrow (2·ζ·ω_n = 2 rad/s). BT's frequency resolution is roughly π / M rad/sample, so the Hann lag window M must be large enough to see the peak. We also pass a dense custom frequency grid so that the peak is not undersampled by the default bin spacing.

In [3]:

w_grid = np.linspace(0.005, np.pi, 512)
result = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts)

w_grid = np.linspace(0.005, np.pi, 512) result = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts)

Plot Bode diagram¶

The magnitude should settle near the DC value 1/k = 0.01 below resonance, peak near ω_n = 10 rad/s (≈ 1.6 Hz), and then roll off at −40 dB/decade. The phase drops through the resonance. The shaded confidence bands come from the estimated noise spectrum propagated through the BT formulas.

In [4]:

fig, axes = plt.subplots(2, 1, figsize=(8, 6))
sid.bode_plot(result)
plt.tight_layout()
plt.show()

fig, axes = plt.subplots(2, 1, figsize=(8, 6)) sid.bode_plot(result) plt.tight_layout() plt.show()

Plot noise spectrum¶

In [5]:

fig, ax = plt.subplots(figsize=(8, 4))
sid.spectrum_plot(result)
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(8, 4)) sid.spectrum_plot(result) plt.tight_layout() plt.show()

Compare different window sizes¶

The Hann window length M sets the BT frequency resolution. For a narrow resonance like this one, a short window (M = 50) smears the peak so badly that the maximum migrates toward DC, while a long window (M = 300) recovers the peak within half a bin of the true ω_n. This is the classic bias/variance trade-off — long windows resolve narrow features but are noisier at smooth parts of the spectrum.

In [6]:

r50  = sid.freq_bt(y, u, window_size=50,  frequencies=w_grid, sample_time=Ts)
r100 = sid.freq_bt(y, u, window_size=100, frequencies=w_grid, sample_time=Ts)
r300 = sid.freq_bt(y, u, window_size=300, frequencies=w_grid, sample_time=Ts)

freq = r300.frequency / Ts   # rad/s

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(freq, 20 * np.log10(np.abs(r50.response)),  'b', label='M = 50')
ax.semilogx(freq, 20 * np.log10(np.abs(r100.response)), 'r', label='M = 100')
ax.semilogx(freq, 20 * np.log10(np.abs(r300.response)), 'g', label='M = 300')
ax.axvline(10.0, color='k', linestyle=':', alpha=0.5, label='true ω_n')
ax.set_xlabel('Frequency (rad/s)')
ax.set_ylabel('Magnitude (dB)')
ax.set_title('Effect of window size on resonance resolution')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

r50 = sid.freq_bt(y, u, window_size=50, frequencies=w_grid, sample_time=Ts) r100 = sid.freq_bt(y, u, window_size=100, frequencies=w_grid, sample_time=Ts) r300 = sid.freq_bt(y, u, window_size=300, frequencies=w_grid, sample_time=Ts) freq = r300.frequency / Ts # rad/s fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogx(freq, 20 * np.log10(np.abs(r50.response)), 'b', label='M = 50') ax.semilogx(freq, 20 * np.log10(np.abs(r100.response)), 'r', label='M = 100') ax.semilogx(freq, 20 * np.log10(np.abs(r300.response)), 'g', label='M = 300') ax.axvline(10.0, color='k', linestyle=':', alpha=0.5, label='true ω_n') ax.set_xlabel('Frequency (rad/s)') ax.set_ylabel('Magnitude (dB)') ax.set_title('Effect of window size on resonance resolution') ax.legend() ax.grid(True) plt.tight_layout() plt.show()

Preprocessing: detrend data before estimation¶

Low-frequency artifacts (sensor drift, DC offsets) bias frequency response estimates near DC. We overlay a linear drift on the output and a DC shift on the input, then show that sid.detrend removes the bias before BT estimation.

In [7]:

y_drift = y + 2e-4 * np.arange(1, N + 1)   # add linear drift (~0.4 m total)
u_drift = u + 5.0                           # add DC offset on input

# Without detrending: drift biases the low-frequency estimate
result_raw = sid.freq_bt(y_drift, u_drift, window_size=200,
                          frequencies=w_grid, sample_time=Ts)

# With detrending
y_dt, _ = sid.detrend(y_drift)
u_dt, _ = sid.detrend(u_drift)
result_dt = sid.freq_bt(y_dt, u_dt, window_size=200,
                         frequencies=w_grid, sample_time=Ts)

print(f'Without detrend: max |G| at low freq = {np.max(np.abs(result_raw.response)):.4f}')
print(f'With detrend:    max |G| at low freq = {np.max(np.abs(result_dt.response)):.4f}')

y_drift = y + 2e-4 * np.arange(1, N + 1) # add linear drift (~0.4 m total) u_drift = u + 5.0 # add DC offset on input # Without detrending: drift biases the low-frequency estimate result_raw = sid.freq_bt(y_drift, u_drift, window_size=200, frequencies=w_grid, sample_time=Ts) # With detrending y_dt, _ = sid.detrend(y_drift) u_dt, _ = sid.detrend(u_drift) result_dt = sid.freq_bt(y_dt, u_dt, window_size=200, frequencies=w_grid, sample_time=Ts) print(f'Without detrend: max |G| at low freq = {np.max(np.abs(result_raw.response)):.4f}') print(f'With detrend: max |G| at low freq = {np.max(np.abs(result_dt.response)):.4f}')

Without detrend: max |G| at low freq = 0.7407
With detrend:    max |G| at low freq = 0.0305

Model validation: residual analysis¶

sid.residual runs two diagnostics on the one-step prediction residuals: a whiteness test (the autocorrelation of the residuals should fall inside a 99% confidence bound) and an independence test (the cross-correlation between residuals and input should also fall inside the bound).

Expected result on this plant: whiteness FAILS. The non-parametric Blackman-Tukey estimator has finite-window bias, and on a narrow resonance (Q = 5 here) that bias leaks into the one-step residuals as a small but systematic auto-correlation that exceeds the 99% confidence bound. This is not a bug in the example — it is the estimator telling you "your model is missing structure that a parametric method could capture". See example_ltv_disc for a parametric COSMIC fit on a related LTV plant where the whiteness test passes cleanly.

In [8]:

resid = sid.residual(result, y, u)

# Print both verdicts; whiteness FAIL is expected on this narrow
# resonance (see the markdown cell above for the explanation).
print('Whiteness test:   ', 'PASS' if resid.whiteness_pass    else 'FAIL (expected for BT on narrow resonances)')
print('Independence test:', 'PASS' if resid.independence_pass else 'FAIL')

resid = sid.residual(result, y, u) # Print both verdicts; whiteness FAIL is expected on this narrow # resonance (see the markdown cell above for the explanation). print('Whiteness test: ', 'PASS' if resid.whiteness_pass else 'FAIL (expected for BT on narrow resonances)') print('Independence test:', 'PASS' if resid.independence_pass else 'FAIL')

Whiteness test:    FAIL (expected for BT on narrow resonances)
Independence test: PASS

Model validation: compare predicted vs measured¶

In [9]:

comp = sid.compare(result, y, u)
print(f'NRMSE fit: {comp.fit[0]:.1f}%')

comp = sid.compare(result, y, u) print(f'NRMSE fit: {comp.fit[0]:.1f}%')

NRMSE fit: 67.3%

Time-series mode (no input)¶

If the input record is unavailable, sid.freq_bt can still estimate the output power spectrum. We re-simulate the SMD under fresh white-force excitation and hand only the position trajectory to BT — the resonance peak should still appear at ω_n ≈ 10 rad/s because the plant colours the white input with its own frequency response.

In [10]:

u_ts = rng.standard_normal(1000)
x_ts = np.zeros((1001, 2))
for step in range(1000):
    x_ts[step + 1] = Ad @ x_ts[step] + Bd[:, 0] * u_ts[step]
y_ts = x_ts[1:, 0]

result_ts = sid.freq_bt(y_ts, None, window_size=200, frequencies=w_grid)

sid.spectrum_plot(result_ts)
plt.title('SDOF output spectrum (time-series mode)')
plt.tight_layout()
plt.show()

u_ts = rng.standard_normal(1000) x_ts = np.zeros((1001, 2)) for step in range(1000): x_ts[step + 1] = Ad @ x_ts[step] + Bd[:, 0] * u_ts[step] y_ts = x_ts[1:, 0] result_ts = sid.freq_bt(y_ts, None, window_size=200, frequencies=w_grid) sid.spectrum_plot(result_ts) plt.title('SDOF output spectrum (time-series mode)') plt.tight_layout() plt.show()