Example: Frequency-Dependent Resolution¶
This example demonstrates sid.freq_btfdr, which uses a different
Hann lag window size at each frequency. This is valuable when a
system has sharp features (like structural resonances) at some
frequencies but is featureless elsewhere — you want fine resolution
only where it matters. Replaces MATLAB's spafdr.
Plant. A 3-mass spring-damper chain
wall--k₁--m₁--k₂--m₂--k₃--m₃ with nonuniform stiffness and
lightly-damped modes:
| Parameter | Value |
|---|---|
m |
[1, 1, 1] kg |
k |
[300, 200, 100] N/m |
c |
[8, 8, 8] N·s/m |
The three normal modes sit at approximately 6.4, 15.1, and
25.1 rad/s — well-separated, all in the bottom decade of the
spectrum. Above the third mode the response rolls off smoothly,
which is exactly the profile where BTFDR shines: fine resolution
around the peaks, coarse resolution in the featureless high
frequencies.
Translated from exampleFreqDepRes.m.
import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd
import sid
%matplotlib inline
Generate test data¶
Force is applied to mass 1 and we measure mass 1's position. We
also keep the true discrete transfer function
G(e^{jω}) = C·(e^{jω}·I − Ad)^{-1}·Bd as a dashed reference
throughout the comparisons below.
rng = np.random.default_rng(2)
# ---- 3-mass SMD chain ----
m = np.array([1.0, 1.0, 1.0])
k = np.array([300.0, 200.0, 100.0])
c = np.array([8.0, 8.0, 8.0])
F = np.array([[1.0], [0.0], [0.0]]) # force on mass 1
Ts = 0.01
N = 6000
Ad, Bd = util_msd(m, k, c, F, Ts)
C_out = np.zeros((1, 6)); C_out[0, 0] = 1.0 # measure position x₁
u = rng.standard_normal(N)
x = np.zeros((N + 1, 6))
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)
def true_tf(w):
'''Exact discrete transfer function at digital frequency w (rad/sample).'''
return np.array([
(C_out @ np.linalg.inv(np.exp(1j * wi) * np.eye(6) - Ad) @ Bd)[0, 0]
for wi in w
])
Fixed-window freq_bt: the resolution-variance trade-off¶
A short Hann window (M = 15) gives a smooth estimate but fails to
distinguish the three modes — they bleed into a single broad blob.
A long window (M = 80) resolves the peaks but the estimate is
noisy in the flat regions above 1 rad/sample.
r_small = sid.freq_bt(y, u, window_size=15, sample_time=Ts)
r_large = sid.freq_bt(y, u, window_size=80, sample_time=Ts)
w = r_small.frequency # rad/sample
G_true = true_tf(w)
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w, 20 * np.log10(np.abs(r_small.response)), 'b', label='BT M = 15')
ax.semilogx(w, 20 * np.log10(np.abs(r_large.response)), 'r', label='BT M = 80')
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('Fixed window: resolution vs variance trade-off')
ax.legend(loc='lower left')
ax.grid(True)
plt.tight_layout()
plt.show()
Scalar resolution with freq_btfdr¶
Passing a scalar resolution = R sets the window size uniformly
from M = ceil(2π / R). Smaller R = finer resolution = larger
window.
result_fdr = sid.freq_btfdr(y, u, resolution=0.2, sample_time=Ts)
sid.bode_plot(result_fdr)
plt.suptitle('freq_btfdr with scalar resolution R = 0.2', y=1.02)
plt.tight_layout()
plt.show()
Per-frequency resolution vector¶
The real power of BTFDR is passing a vector resolution: fine
at frequencies where the plant has sharp features, coarse where it
is smooth. For our 3-mass plant, the modes are in the bottom
decade, so we ramp R from 0.1 (fine) at low frequencies to
1.5 (coarse) at high frequencies.
nf = len(result_fdr.frequency)
R_vec = np.linspace(0.1, 1.5, nf) # fine at low freq, coarse at high freq
result_vec = sid.freq_btfdr(y, u, resolution=R_vec, sample_time=Ts)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7))
ax1.plot(result_vec.frequency, result_vec.window_size, 'b')
ax1.set_xlabel('Frequency (rad/sample)')
ax1.set_ylabel('Window size M')
ax1.set_title('Per-frequency window size')
ax1.grid(True)
ax2.semilogx(result_vec.frequency,
20 * np.log10(np.abs(result_vec.response)),
'b', label='BTFDR (variable R)')
ax2.semilogx(w, 20 * np.log10(np.abs(G_true)),
'k--', lw=1.5, label='True')
ax2.set_xlabel('Frequency (rad/sample)')
ax2.set_ylabel('Magnitude (dB)')
ax2.set_title('BTFDR with per-frequency resolution')
ax2.legend(loc='lower left')
ax2.grid(True)
plt.tight_layout()
plt.show()
Compare BT vs BTFDR side by side¶
BTFDR adapts: it uses a big window near the modes to resolve the peaks, and a small window at high frequencies where the response is featureless, reducing variance there.
r_bt = sid.freq_bt(y, u, window_size=30, sample_time=Ts)
r_fdr = sid.freq_btfdr(y, u, resolution=0.3, sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(r_bt.frequency, 20 * np.log10(np.abs(r_bt.response)),
'b', label='BT (M = 30)')
ax.semilogx(r_fdr.frequency, 20 * np.log10(np.abs(r_fdr.response)),
'r', label='BTFDR (R = 0.3)')
ax.semilogx(r_bt.frequency, 20 * np.log10(np.abs(true_tf(r_bt.frequency))),
'k--', lw=1.5, label='True')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_ylabel('Magnitude (dB)')
ax.set_title('Blackman-Tukey: fixed vs frequency-dependent resolution')
ax.legend(loc='lower left')
ax.grid(True)
plt.tight_layout()
plt.show()
