Example: Coherence Analysis and Signal Quality Assessment¶
This example uses the squared coherence γ²(ω) to assess the
quality of a frequency response estimate. Coherence near 1 means
the output at frequency ω is explained well by the input; values
near 0 mean noise or unmodeled disturbances dominate. Coherence is
only available for SISO estimates from sid.freq_bt or
sid.freq_btfdr (not sid.freq_etfe).
Plant. A 2-mass spring-damper chain with a structural disturbance entering at the second mass:
| Parameter | Value |
|---|---|
m |
[1, 1] kg |
k |
[100, 80] N/m |
c |
[2, 2] N·s/m |
The operator applies a commanded force u(t) at mass 1 and
measures mass 1's position. A secondary, uncommanded force d(t)
enters at mass 2 — think of an unmodeled actuator, a nearby
machine, or a sloshing fluid. d(t) is colored (low-frequency
dominant, AR(1) pole at 0.9), so the coherence between u and
y should be worst at low frequencies where d has the most
power. Normal modes are at ω ≈ 6.0 and 15.0 rad/s.
Translated from exampleCoherence.m.
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lfilter
from util_msd import util_msd
import sid
%matplotlib inline
Generate test data¶
The F matrix has two columns: column 0 injects the commanded
force u at mass 1; column 1 injects the disturbance d at
mass 2. util_msd returns a single Bd of shape (4, 2) that
handles both channels simultaneously.
rng = np.random.default_rng(3)
# ---- 2-mass SMD with 2 input channels (command + disturbance) ----
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],
[0.0, 1.0]]) # col 0: u at mass 1; col 1: d at mass 2
Ts = 0.01
N = 4000
Ad, Bd = util_msd(m, k, c, F, Ts)
C_out = np.zeros((1, 4)); C_out[0, 0] = 1.0 # measure mass 1 position
# ---- Excitation and disturbance ----
u = 10.0 * rng.standard_normal(N) # commanded force
e = rng.standard_normal(N)
d = 0.5 * lfilter([1], [1, -0.9], e) # DC-heavy colored disturbance
# ---- Simulate ----
x = np.zeros((N + 1, 4))
for step in range(N):
x[step + 1] = Ad @ x[step] + Bd[:, 0] * u[step] + Bd[:, 1] * d[step]
y = x[1:, 0]
Estimate with freq_bt¶
A custom frequency grid is supplied so the narrow modes of the lightly-damped chain are not undersampled.
w_grid = np.linspace(0.01, np.pi, 512)
result = sid.freq_bt(y, u, window_size=200, frequencies=w_grid, sample_time=Ts)
Plot Bode magnitude and coherence together¶
Coherence drops wherever unmodeled content (disturbance or noise floor) masks the cause-effect relationship: near the low-frequency disturbance band, at the antiresonance between the two modes (where the plant gain itself is tiny), and in the high-frequency tail (where plant roll-off lets noise take over).
w = result.frequency
coh = result.coherence.ravel()
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
ax1.semilogx(w, 20 * np.log10(np.abs(result.response.ravel())), 'b')
ax1.set_ylabel('Magnitude (dB)')
ax1.set_title('Frequency response estimate')
ax1.grid(True)
ax2.semilogx(w, coh, 'b')
ax2.axhline(0.9, color='g', ls='--', label=r'$\gamma^2 = 0.9$')
ax2.axhline(0.5, color='r', ls='--', label=r'$\gamma^2 = 0.5$')
ax2.set_ylabel(r'Coherence $\gamma^2$')
ax2.set_xlabel('Frequency (rad/sample)')
ax2.set_title('Squared coherence')
ax2.legend(loc='lower left')
ax2.set_ylim([0, 1])
ax2.grid(True)
plt.tight_layout()
plt.show()
Confidence bands reflect coherence¶
Where coherence is high the uncertainty is low (narrow bands).
The confidence option controls the number of standard
deviations shown around the estimate.
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
sid.bode_plot(result, confidence=2, ax=(axes[0, 0], axes[1, 0]))
axes[0, 0].set_title(r'$2\sigma$ confidence bands')
sid.bode_plot(result, confidence=3, ax=(axes[0, 1], axes[1, 1]))
axes[0, 1].set_title(r'$3\sigma$ confidence bands')
plt.tight_layout()
plt.show()
High-disturbance vs low-disturbance comparison¶
Scaling the colored disturbance up drags coherence down everywhere the disturbance has spectral power. Scaling it down recovers coherence near the plant modes where the signal dominates.
rng2 = np.random.default_rng(3)
e2 = rng2.standard_normal(N)
d_low = 0.1 * lfilter([1], [1, -0.9], e2)
d_high = 2.0 * lfilter([1], [1, -0.9], e2)
def simulate(d_signal):
xs = np.zeros((N + 1, 4))
for step in range(N):
xs[step + 1] = Ad @ xs[step] + Bd[:, 0] * u[step] + Bd[:, 1] * d_signal[step]
return xs[1:, 0]
y_low = simulate(d_low)
y_high = simulate(d_high)
r_low = sid.freq_bt(y_low, u, window_size=200, frequencies=w_grid, sample_time=Ts)
r_high = sid.freq_bt(y_high, u, window_size=200, frequencies=w_grid, sample_time=Ts)
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(w, r_low.coherence.ravel(), 'b', label=r'Low disturbance (0.1 $\times$)')
ax.semilogx(w, r_high.coherence.ravel(), 'r', label=r'High disturbance (2.0 $\times$)')
ax.set_ylabel(r'Coherence $\gamma^2$')
ax.set_xlabel('Frequency (rad/sample)')
ax.set_title('Effect of disturbance level on coherence')
ax.legend(loc='lower left')
ax.set_ylim([0, 1])
ax.grid(True)
plt.tight_layout()
plt.show()
Note: ETFE does not provide coherence¶
sid.freq_etfe returns coherence = None because there is no
windowed cross-spectral estimate in the FFT-ratio approach.
result_etfe = sid.freq_etfe(y, u)
print(f'ETFE coherence is None: {result_etfe.coherence is None}')
ETFE coherence is None: True