Example: Empirical Transfer Function Estimate (ETFE)¶
This example demonstrates sid.freq_etfe, which estimates the
frequency response as the ratio of output and input DFTs. It has
maximum frequency resolution but high variance at every bin;
optional smoothing averages adjacent bins to trade resolution for
variance.
Plant. The same 1-DoF spring-mass-damper as example_siso:
m = 1 kg, k = 100 N/m, c = 2 N·s/m — natural frequency
ω_n = 10 rad/s, damping ratio ζ = 0.1. Input is a force, output
is the mass position x₁. Using the same plant across both
notebooks lets us compare how BT and ETFE handle the same physical
problem.
Translated from exampleETFE.m.
import numpy as np
import matplotlib.pyplot as plt
from util_msd import util_msd
import sid
%matplotlib inline
Generate test data¶
Build the SDOF state-space model with util_msd, simulate under
white-force excitation, and add sensor noise to the measured
position. We also compute the true discrete transfer function
G(e^{jω}) = C·(e^{jω}·I − Ad)^{-1}·Bd as a reference for the
smoothing comparison below.
rng = np.random.default_rng(1)
# ---- Physical plant: 1-DoF SMD (ω_n = 10 rad/s, ζ = 0.1) ----
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]]) # measure position
# ---- Simulate the plant ----
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_clean = x[1:, 0]
y = y_clean + 2e-4 * rng.standard_normal(N)
Basic ETFE (no smoothing)¶
The raw ETFE is just Y(ω) / U(ω) — maximum resolution, extremely
noisy. response_std is NaN because ETFE has no asymptotic
uncertainty formula.
result = sid.freq_etfe(y, u, sample_time=Ts)
sid.bode_plot(result)
plt.suptitle('ETFE - raw (no smoothing)', y=1.02)
plt.tight_layout()
plt.show()
Effect of smoothing¶
Smoothing averages S adjacent frequency bins (S must be odd).
Larger S suppresses the periodogram noise but smears the
resonance peak. We overlay the true discrete transfer function as a
dashed reference.
r1 = sid.freq_etfe(y, u, smoothing=1, sample_time=Ts) # raw
r11 = sid.freq_etfe(y, u, smoothing=11, sample_time=Ts) # moderate
r21 = sid.freq_etfe(y, u, smoothing=21, sample_time=Ts) # heavy
w = r1.frequency # rad/sample
# True TF at the ETFE bin grid
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(r1.response)), color='0.7', label='S = 1 (raw)')
ax.semilogx(w, 20 * np.log10(np.abs(r11.response)), 'b', label='S = 11')
ax.semilogx(w, 20 * np.log10(np.abs(r21.response)), 'r', label='S = 21')
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('ETFE smoothing comparison')
ax.legend(loc='lower left')
ax.grid(True)
plt.tight_layout()
plt.show()
Known FIR system: pure delay¶
For y(t) = u(t − 1) the true transfer function is G(z) = z^{-1}:
unit magnitude and phase -ω at every frequency. ETFE recovers this
exactly when there is no noise. This is a reality check on the
estimator independent of the SMD plant above.
rng2 = np.random.default_rng(2)
N_fir = 1024
u_fir = rng2.standard_normal(N_fir)
y_fir = np.concatenate(([0.0], u_fir[:-1])) # one-sample delay
result_fir = sid.freq_etfe(y_fir, u_fir)
w_fir = result_fir.frequency
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6))
ax1.plot(w_fir, np.abs(result_fir.response), 'b')
ax1.set_ylabel('|G|')
ax1.set_title('ETFE of pure delay: |G| should be 1')
ax1.grid(True)
ax2.plot(w_fir, np.angle(result_fir.response), 'b', label='ETFE')
ax2.plot(w_fir, -w_fir, 'k--', label=r'True phase = $-\omega$')
ax2.set_ylabel('Phase (rad)')
ax2.set_xlabel('Frequency (rad/sample)')
ax2.legend(loc='lower left')
ax2.grid(True)
plt.tight_layout()
plt.show()
Time-series mode: periodogram¶
With no input signal, sid.freq_etfe computes the periodogram of
the output. We re-simulate the same SDOF under fresh white-force
excitation and pass only the position record. The resonance still
appears as a peak near ω_n because the plant colours the input's
flat spectrum with its own frequency response.
rng3 = np.random.default_rng(3)
N_ts = 500
u_ts = rng3.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]
result_ts = sid.freq_etfe(y_ts, None)
sid.spectrum_plot(result_ts)
plt.title('SDOF output periodogram')
plt.tight_layout()
plt.show()
Custom frequency grid and Hz display¶
For narrow resonances the uniform DFT grid can undersample the
peak. Supply a custom log-spaced grid and plot in Hz using
bode_plot(..., frequency_unit='Hz').
rng4 = np.random.default_rng(4)
# Log-spaced grid in rad/sample, spanning 0.005 to π
w_log = np.logspace(np.log10(0.005), np.log10(np.pi), 200)
N_hz = 2048
u_hz = rng4.standard_normal(N_hz)
x_hz = np.zeros((N_hz + 1, 2))
for step in range(N_hz):
x_hz[step + 1] = Ad @ x_hz[step] + Bd[:, 0] * u_hz[step]
y_hz = x_hz[1:, 0] + 2e-4 * rng4.standard_normal(N_hz)
result_hz = sid.freq_etfe(y_hz, u_hz, smoothing=11,
frequencies=w_log, sample_time=Ts)
sid.bode_plot(result_hz, frequency_unit='Hz')
plt.suptitle('ETFE with log frequency grid (Hz)', y=1.02)
plt.tight_layout()
plt.show()
