Example: Comparing Frequency Response Estimation Methods¶

This example compares the three frequency-domain estimators in sid on the same physical plant used by example_siso and example_etfe — a 1-DoF spring-mass-damper with a lightly-damped resonance (ω_n = 10 rad/s, ζ = 0.1, Q = 5):

• sid.freq_bt — Blackman-Tukey (fixed window, replaces spa)
• sid.freq_btfdr — Blackman-Tukey with frequency-dependent resolution (replaces spafdr)
• sid.freq_etfe — Empirical Transfer Function Estimate (replaces etfe)

A preview of the result: on a narrow resonance, the raw ETFE (full frequency resolution, high variance) can actually beat the smoothed Blackman-Tukey estimate on NRMSE fit. Smoothing the ETFE or shrinking the BT window erases the peak and costs fit. This cuts against the common wisdom that BT is always the safe default — for lightly-damped structures, maximum resolution wins.

Translated from exampleMethodComparison.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¶

Reuse the 1-DoF SMD plant from example_siso. A 2048-sample record with unit-variance white-force excitation and 0.2 mm sensor noise.

In [2]:

rng = np.random.default_rng(4)

# ---- 1-DoF SMD plant (same as example_siso) ----
m  = np.array([1.0])
k  = np.array([100.0])
c  = np.array([2.0])
F  = np.array([[1.0]])
Ts = 0.01
N  = 2048

Ad, Bd = util_msd(m, k, c, F, Ts)
C_out = np.array([[1.0, 0.0]])

u = rng.standard_normal(N)
x = np.zeros((N + 1, 2))
for step in range(N):
    x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step]
y = x[1:, 0] + 2e-4 * rng.standard_normal(N)

rng = np.random.default_rng(4) # ---- 1-DoF SMD plant (same as example_siso) ---- m = np.array([1.0]) k = np.array([100.0]) c = np.array([2.0]) F = np.array([[1.0]]) Ts = 0.01 N = 2048 Ad, Bd = util_msd(m, k, c, F, Ts) C_out = np.array([[1.0, 0.0]]) u = rng.standard_normal(N) x = np.zeros((N + 1, 2)) for step in range(N): x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step] y = x[1:, 0] + 2e-4 * rng.standard_normal(N)

Estimate with all three methods¶

All three calls share a custom frequency grid (linspace(0.005, π, 512)) so they compare on identical bins.

In [3]:

w_grid = np.linspace(0.005, np.pi, 512)

r_bt     = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts)
r_etfe   = sid.freq_etfe(y, u, frequencies=w_grid, sample_time=Ts)
r_etfe_s = sid.freq_etfe(y, u, smoothing=15, frequencies=w_grid, sample_time=Ts)
r_fdr    = sid.freq_btfdr(y, u, resolution=0.3, frequencies=w_grid, sample_time=Ts)

w_grid = np.linspace(0.005, np.pi, 512) r_bt = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts) r_etfe = sid.freq_etfe(y, u, frequencies=w_grid, sample_time=Ts) r_etfe_s = sid.freq_etfe(y, u, smoothing=15, frequencies=w_grid, sample_time=Ts) r_fdr = sid.freq_btfdr(y, u, resolution=0.3, frequencies=w_grid, sample_time=Ts)

Compare Bode magnitude plots¶

The dashed reference is the exact discrete transfer function C·(e^{jω}·I − Ad)^{-1}·Bd.

In [4]:

w = r_bt.frequency
G_true = np.array([
    (C_out @ np.linalg.inv(np.exp(1j * wi) * np.eye(2) - Ad) @ Bd)[0, 0]
    for wi in w
])

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w, 20 * np.log10(np.abs(r_etfe.response)),
            color='0.8', label='ETFE (raw)')
ax.semilogx(w, 20 * np.log10(np.abs(r_etfe_s.response)),
            'g', label='ETFE (S = 15)')
ax.semilogx(w, 20 * np.log10(np.abs(r_bt.response)),
            'b', label='BT (M = 200)')
ax.semilogx(w, 20 * np.log10(np.abs(r_fdr.response)),
            'r', label='BTFDR (R = 0.3)')
ax.semilogx(w, 20 * np.log10(np.abs(G_true)),
            'k--', lw=1.5, label='True')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel('Magnitude (dB)')
ax.set_title('Method comparison: Bode magnitude')
ax.legend(loc='lower left')
ax.grid(True)
plt.tight_layout()
plt.show()

w = r_bt.frequency G_true = np.array([ (C_out @ np.linalg.inv(np.exp(1j * wi) * np.eye(2) - Ad) @ Bd)[0, 0] for wi in w ]) fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogx(w, 20 * np.log10(np.abs(r_etfe.response)), color='0.8', label='ETFE (raw)') ax.semilogx(w, 20 * np.log10(np.abs(r_etfe_s.response)), 'g', label='ETFE (S = 15)') ax.semilogx(w, 20 * np.log10(np.abs(r_bt.response)), 'b', label='BT (M = 200)') ax.semilogx(w, 20 * np.log10(np.abs(r_fdr.response)), 'r', label='BTFDR (R = 0.3)') ax.semilogx(w, 20 * np.log10(np.abs(G_true)), 'k--', lw=1.5, label='True') ax.set_xlabel('Frequency (rad/sample)') ax.set_ylabel('Magnitude (dB)') ax.set_title('Method comparison: Bode magnitude') ax.legend(loc='lower left') ax.grid(True) plt.tight_layout() plt.show()

Compare noise spectra¶

BT and BTFDR compute the noise spectrum from covariance estimates, ETFE from residuals. They should all trace out roughly the same noise-floor shape.

In [5]:

# Clamp the noise spectrum with a small floor before log10: per
# SPEC.md 2.7 the noise spectrum is non-negative and may hit exact
# zero at some frequencies (negative estimator variance clamped to
# 0), which would otherwise produce log10(0) warnings.
floor = 1e-30

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w, 10 * np.log10(np.maximum(np.abs(r_bt.noise_spectrum), floor)),
            'b', label='BT')
ax.semilogx(w, 10 * np.log10(np.maximum(np.abs(r_etfe_s.noise_spectrum), floor)),
            'g', label='ETFE (S = 15)')
ax.semilogx(w, 10 * np.log10(np.maximum(np.abs(r_fdr.noise_spectrum), floor)),
            'r', label='BTFDR')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel('Noise spectrum (dB)')
ax.set_title('Noise spectrum comparison')
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

# Clamp the noise spectrum with a small floor before log10: per # SPEC.md 2.7 the noise spectrum is non-negative and may hit exact # zero at some frequencies (negative estimator variance clamped to # 0), which would otherwise produce log10(0) warnings. floor = 1e-30 fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogx(w, 10 * np.log10(np.maximum(np.abs(r_bt.noise_spectrum), floor)), 'b', label='BT') ax.semilogx(w, 10 * np.log10(np.maximum(np.abs(r_etfe_s.noise_spectrum), floor)), 'g', label='ETFE (S = 15)') ax.semilogx(w, 10 * np.log10(np.maximum(np.abs(r_fdr.noise_spectrum), floor)), 'r', label='BTFDR') ax.set_xlabel('Frequency (rad/sample)') ax.set_ylabel('Noise spectrum (dB)') ax.set_title('Noise spectrum comparison') ax.legend() ax.grid(True) plt.tight_layout() plt.show()

Custom logarithmic frequency grid¶

A log-spaced grid covers the low-frequency region more densely without wasting bins on the featureless high-frequency tail.

In [6]:

w_log = np.logspace(np.log10(0.005), np.log10(np.pi), 200)

r_bt_log   = sid.freq_bt(y, u, window_size=200, frequencies=w_log, sample_time=Ts)
r_etfe_log = sid.freq_etfe(y, u, smoothing=15, frequencies=w_log, sample_time=Ts)
r_fdr_log  = sid.freq_btfdr(y, u, resolution=0.3, frequencies=w_log, sample_time=Ts)

G_true_log = np.array([
    (C_out @ np.linalg.inv(np.exp(1j * wi) * np.eye(2) - Ad) @ Bd)[0, 0]
    for wi in w_log
])

fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w_log, 20 * np.log10(np.abs(r_bt_log.response)),
            'b', label='BT')
ax.semilogx(w_log, 20 * np.log10(np.abs(r_etfe_log.response)),
            'g', label='ETFE')
ax.semilogx(w_log, 20 * np.log10(np.abs(r_fdr_log.response)),
            'r', label='BTFDR')
ax.semilogx(w_log, 20 * np.log10(np.abs(G_true_log)),
            'k--', lw=1.5, label='True')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel('Magnitude (dB)')
ax.set_title('Log frequency grid (200 points)')
ax.legend(loc='lower left')
ax.grid(True)
plt.tight_layout()
plt.show()

w_log = np.logspace(np.log10(0.005), np.log10(np.pi), 200) r_bt_log = sid.freq_bt(y, u, window_size=200, frequencies=w_log, sample_time=Ts) r_etfe_log = sid.freq_etfe(y, u, smoothing=15, frequencies=w_log, sample_time=Ts) r_fdr_log = sid.freq_btfdr(y, u, resolution=0.3, frequencies=w_log, sample_time=Ts) G_true_log = np.array([ (C_out @ np.linalg.inv(np.exp(1j * wi) * np.eye(2) - Ad) @ Bd)[0, 0] for wi in w_log ]) fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogx(w_log, 20 * np.log10(np.abs(r_bt_log.response)), 'b', label='BT') ax.semilogx(w_log, 20 * np.log10(np.abs(r_etfe_log.response)), 'g', label='ETFE') ax.semilogx(w_log, 20 * np.log10(np.abs(r_fdr_log.response)), 'r', label='BTFDR') ax.semilogx(w_log, 20 * np.log10(np.abs(G_true_log)), 'k--', lw=1.5, label='True') ax.set_xlabel('Frequency (rad/sample)') ax.set_ylabel('Magnitude (dB)') ax.set_title('Log frequency grid (200 points)') ax.legend(loc='lower left') ax.grid(True) plt.tight_layout() plt.show()

Time-series comparison: periodogram vs smoothed spectrum¶

With no input (u = None), all methods estimate the output power spectrum of the same SMD driven by fresh white force noise. The raw ETFE periodogram is noisy; BT is much smoother.

In [7]:

rng5 = np.random.default_rng(5)

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

r_ts_bt   = sid.freq_bt(y_ts, None, window_size=100)
r_ts_etfe = sid.freq_etfe(y_ts, None)

w_ts = r_ts_bt.frequency
floor = 1e-30
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w_ts, 10 * np.log10(np.maximum(np.abs(r_ts_etfe.noise_spectrum), floor)),
            color='0.7', label='ETFE (periodogram)')
ax.semilogx(w_ts, 10 * np.log10(np.maximum(np.abs(r_ts_bt.noise_spectrum), floor)),
            'b', label='BT (M = 100)')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel('Power spectrum (dB)')
ax.set_title('Time-series: periodogram vs Blackman-Tukey')
ax.legend(loc='lower left')
ax.grid(True)
plt.tight_layout()
plt.show()

rng5 = np.random.default_rng(5) N_ts = 1000 u_ts = rng5.standard_normal(N_ts) x_ts = np.zeros((N_ts + 1, 2)) for step in range(N_ts): x_ts[step + 1] = Ad @ x_ts[step] + Bd[:, 0] * u_ts[step] y_ts = x_ts[1:, 0] r_ts_bt = sid.freq_bt(y_ts, None, window_size=100) r_ts_etfe = sid.freq_etfe(y_ts, None) w_ts = r_ts_bt.frequency floor = 1e-30 fig, ax = plt.subplots(figsize=(8, 5)) ax.semilogx(w_ts, 10 * np.log10(np.maximum(np.abs(r_ts_etfe.noise_spectrum), floor)), color='0.7', label='ETFE (periodogram)') ax.semilogx(w_ts, 10 * np.log10(np.maximum(np.abs(r_ts_bt.noise_spectrum), floor)), 'b', label='BT (M = 100)') ax.set_xlabel('Frequency (rad/sample)') ax.set_ylabel('Power spectrum (dB)') ax.set_title('Time-series: periodogram vs Blackman-Tukey') ax.legend(loc='lower left') ax.grid(True) plt.tight_layout() plt.show()

Model output comparison using sid.compare¶

For each estimator, simulate the forward model back through the measured input and report the NRMSE fit to the measured output. On this lightly-damped plant, the raw ETFE wins by a large margin — it preserves the full resolution of the resonance that BT/BTFDR/smoothed ETFE all smear.

In [8]:

comp_bt     = sid.compare(r_bt,     y, u)
comp_fdr    = sid.compare(r_fdr,    y, u)
comp_etfe   = sid.compare(r_etfe,   y, u)
comp_etfe_s = sid.compare(r_etfe_s, y, u)

print('--- Prediction fit (NRMSE %) ---')
print(f'  freq_bt    (M=200):  {comp_bt.fit[0]:5.1f}%')
print(f'  freq_btfdr (R=0.3):  {comp_fdr.fit[0]:5.1f}%')
print(f'  freq_etfe  (raw):    {comp_etfe.fit[0]:5.1f}%')
print(f'  freq_etfe  (S=15):   {comp_etfe_s.fit[0]:5.1f}%')

comp_bt = sid.compare(r_bt, y, u) comp_fdr = sid.compare(r_fdr, y, u) comp_etfe = sid.compare(r_etfe, y, u) comp_etfe_s = sid.compare(r_etfe_s, y, u) print('--- Prediction fit (NRMSE %) ---') print(f' freq_bt (M=200): {comp_bt.fit[0]:5.1f}%') print(f' freq_btfdr (R=0.3): {comp_fdr.fit[0]:5.1f}%') print(f' freq_etfe (raw): {comp_etfe.fit[0]:5.1f}%') print(f' freq_etfe (S=15): {comp_etfe_s.fit[0]:5.1f}%')

--- Prediction fit (NRMSE %) ---
  freq_bt    (M=200):   63.4%
  freq_btfdr (R=0.3):   10.1%
  freq_etfe  (raw):     91.2%
  freq_etfe  (S=15):    32.5%

Summary of method trade-offs¶

Method	Window size	Uncertainty	Coherence	Best for
freq_bt	Fixed M	Yes	Yes	Broad features, smooth plants
freq_btfdr	Per-frequency M_k	Yes	Yes	Plants with features at known frequency bands
freq_etfe	N (full DFT)	No	No	Narrow resonances; preliminary exploration

freq_bt remains the right default when you have no prior idea where the plant features sit and want uncertainty bands on every estimate. freq_etfe raw is the right tool for sharp resonances where resolution matters more than smoothness. freq_btfdr is worth the extra configuration when you know the rough band structure of your plant.