# Method 7 — Output-only Closed-loop Resonance / Sensitivity-peak NCF Margin

*Conclusions, rationale, and a detailed implementation guide.*

---

## 0. TL;DR

When the reference input does **not** excite the margin-binding frequency
(~35 Hz here), the open-loop FRF `H_dc_dm = dm/dc` is **uninformative there** —
it is pinned at `1/|C|` by the feedback law, no matter how close to instability
the loop is. But the **closed-loop resonance still grows in the response `dm`
alone**, because near instability the loop amplifies internal disturbance at the
resonance. Method 7 reads that growth:

```
Spk(t) = sqrt( peak PSD of dm in [25,40] Hz / floor PSD of dm in [42,58] Hz )   ~  ||S||_inf
b(t)   = b_nom * Spk_ref / Spk(t)        % calibrated to the model's nominal margin at the start
```

It is **input-agnostic**: it catches the instability on the low-frequency-input
file (1749, detect 219 s; b → 0.07 at 256–270 s — the *only* method that does)
**and** agrees with the FRF methods on the broadband file (1731, detect 717 s vs
their 728–734 s). Implementation: `common/osi_method7_resonance.m`, runner
`trial_7/run_trial_7.m`, parameters `cfg.m7_*` in `common/osi_config.m`.

---

## 1. What problem made me propose Method 7

Methods 1–6 estimate the loop `L = C·H_dc_dm` from a windowed FRF and take the
NCF margin `b = min_ω |1+L| / sqrt(2(1+|L|²))`. They work beautifully on the
broadband record (1731): the margin tracks a clean collapse and trips ~715–734 s.

On the **low-frequency-input** record (1749) they failed. The instability is the
**same physical mechanism** — the 35 Hz closed-loop resonance growing — but the
input `dd` has a median frequency of ~2.8 Hz and essentially **no energy above
~6 Hz**. My first (wrong) reaction was a power gate that discarded the 35 Hz band
and a conclusion that the mode was "unobservable." It is not. The user correctly
pointed out the 35 Hz resonance is plainly visible and growing in `dm`.

## 2. The key discovery — a closed-loop *identifiability tautology*

The decisive question: *can `H_dc_dm = dm/dc` see the growing 35 Hz resonance?*

**Theory.** The controller is a hard, known constraint: `dc = C·(dd − dm)`. At a
frequency where `dd` has no energy, `dc ≈ −C·dm`, so

```
H_dc_dm(35) = dm/dc ≈ dm / (−C·dm) = −1/C(35)
=>  L = C·H ≈ −1   (a constant)   — independent of the true plant or margin.
```

`dc` and `dm` don’t move independently — they **ring together**, locked by the
controller. The growth near instability is in their **common amplitude**, not in
their **ratio**.

**Empirical proof** (`results/diagnostic_dmdc_vs_amplitude.png`), at 35 Hz:

| | `|dm/dc|` (the FRF) | `dm` amplitude | coherence(dc,dm) |
|---|---|---|---|
| **1749** (low-freq), stable → unstable | **3.10 → 3.10** (= `1/|C|`, flat; *not* the true `|H0|=0.50`) | 0.056 → **3.27** (×60) | **1.00** (perfectly locked) |
| **1731** (broadband), stable → unstable | 0.5 → **2.0** (grows; independent 35 Hz excitation breaks the lock) | 0.12 → 1.0 | 0.4–1.0 |

So on 1749 the FRF ratio carries **zero** information about the instability, while
the response amplitude carries all of it. This is *closed-loop identifiability*,
not a gate I can loosen: you cannot separate plant from controller at a frequency
the input never probes.

## 3. Why the resonance amplitude *is* the margin (the principle)

Write the closed-loop sensitivity `S = 1/(1+L)`. With broadband internal
disturbance `d` (process/sensor noise, nonlinearity) and no reference at 35 Hz,
the response there is `dm ≈ S·d`. Hence the `dm` power spectrum is

```
PSD_dm(ω) ≈ |S(ω)|² · PSD_d(ω)
```

* At the **resonance** `ω_res ≈ 35 Hz`, `|S|` peaks → `peak PSD ∝ |S(ω_res)|²`.
* At an **off-resonance floor** `ω_floor` (here 42–58 Hz, away from the 60 Hz
  line), the loop gain `|L|` is small so `|S| ≈ 1` → `floor PSD ∝ PSD_d`.

Assuming the disturbance is roughly flat across `[ω_res, ω_floor]`, the ratio
removes the unknown disturbance level:

```
Spk = sqrt( peak / floor )  ≈  |S(ω_res)|  =  ||S||_inf  =  1 / min_ω |1+L| .
```

The loop-only NCF near the binding resonance is
`b = min|1+L| / sqrt(2(1+|L|²)) = k · min|1+L| = k / ||S||_inf`, i.e.
**`b ∝ 1/Spk`**. The constant `k` (≈ ½ when `|L(ω_res)| ≈ 1`, the resonance is
near gain crossover) is fixed by the model **calibrated at the start**, where the
margin equals the known nominal `b_nom` and `Spk = Spk_ref`:

```
b(t) = b_nom · Spk_ref / Spk(t).
```

That is Method 7. It estimates the margin from the closed-loop *effect*
(resonance amplification) instead of the *cause* (the plant FRF), which is the
only observable quantity when the input does not excite the binding mode.

## 4. What I learned (the conclusions)

1. **Stability info is in the loop’s closed-loop response, not necessarily in the
   plant FRF.** FRF identification needs persistent excitation at the binding
   frequency; the closed-loop resonance peak does not.
2. **Coherence is *not* a reliability check in closed loop.** On 1749 it reads a
   perfect 1.00 *because* `dc` and `dm` are locked — exactly the case where the
   FRF is meaningless. Use **input power** / external-excitation checks, not
   coherence, to decide whether an FRF point is identifiable.
3. **The two records have the same stability trend**; the only difference is
   *observability* of the binding mode through the chosen signals.
4. **The "model correct at the beginning" assumption is what makes an absolute
   number possible** from a relative measurement: it pins the proportionality
   constant between resonance amplification and NCF margin.
5. **Method 7 is the input-agnostic safety net.** It catches both records, and on
   the broadband one it independently reproduces the FRF methods’ detection — a
   cross-validation that it reads the same physics.

## 5. Scope, assumptions, and limitations

**Use it when** the binding instability shows as a growing resonance in the
measured output (true here, and typical of RTHS actuator/specimen modes), and
especially when the excitation may not cover the binding frequency.

**Assumptions:**
- A single dominant binding resonance in a known band (here 25–40 Hz).
- The internal disturbance is broadband and roughly stationary across the peak
  and floor bands (so the floor normalises it out).
- A clean off-resonance floor band exists (avoid line frequencies — we skip 60 Hz
  by using 42–58 Hz).
- The model is correct at the start (for the calibration).

**Limitations:**
- It returns a **relative / start-calibrated** margin, not an absolute one; the
  constant `k` is empirical.
- It needs a **longer buffer (4 s)** than the FRF methods (1.5 s) for a stable
  spectral peak — coarser time resolution.
- After *full* divergence/saturation the peak/floor ratio can saturate (the whole
  spectrum is the oscillation); it is an approach-to-instability indicator, best
  paired with a raw output-energy trip for the saturated regime.

## 6. The algorithm, step by step

For each update time `t_k` (every `m7_hop = 0.5 s`):

1. **Buffer.** Take the most recent `m7_buf = 4 s` of `dm` (the response only),
   detrended.
2. **Spectrum.** `pwelch(dm, hann(seg), 50% overlap, nfft, fs)` with `seg ≈ 2 s`
   sub-segments → a few averages, ~0.5 Hz resolution.
3. **Resonance amplification.**
   `Spk = sqrt( max PSD over [25,40] Hz / median PSD over [42,58] Hz )`.
4. **Smooth.** `Spk ← medfilt1(Spk, 9)` over the time series (robust to spikes).
5. **Calibrate once** at the start: `Spk_ref = median(Spk over first 20 s)`.
6. **Map to margin.** `b = min( b_nom · Spk_ref / Spk , cap=0.70 )`, where
   `b_nom` = the *model’s* loop-only NCF over the binding band `[10,58] Hz`
   (≈0.588), computed from `C·H0d`.
7. **Detect.** Alarm when `b < 0.40` for `m7_detect_consec = 7` consecutive
   updates (3.5 s sustained).

## 7. Parameters and how to choose them (`cfg.m7_*`)

| Parameter | Default | Meaning / how to set |
|---|---|---|
| `m7_buf` | 4 s | Analysis buffer. Long enough for a stable resonance peak (≥ ~30 cycles of 35 Hz); shorter ⇒ noisier, longer ⇒ laggier. |
| `m7_hop` | 0.5 s | Update period (time resolution). |
| `m7_peak_band` | `[25 40]` Hz | Search band for the binding resonance. Set from a baseline FFT of `dm`; widen if the resonance shifts (it crept to ~37–41 Hz near divergence). |
| `m7_floor_band` | `[42 58]` Hz | Off-resonance noise floor. Must avoid line frequencies (we skip 60 Hz). |
| `m7_binding_band` | `[10 58]` Hz | Band whose **model** NCF gives `b_nom` for the calibration. |
| `m7_calib_s` | 20 s | "Model-correct-at-start" calibration window; keep within the healthy phase. |
| `m7_medfilt` | 9 | Median taps on `Spk` (≈ 4.5 s smoothing on the 0.5 s grid). |
| `m7_cap` | 0.70 | Upper clamp on the reported margin (quiet windows give `Spk < Spk_ref`). |
| `m7_detect_consec` | 7 | Sustained samples below threshold to declare detection (rejects transients). |

## 8. Reference implementation (`common/osi_method7_resonance.m`)

```matlab
function out = osi_method7_resonance(ctx, cfg)
    fs = ctx.fs;
    % nominal binding NCF from the MODEL (calibration target b_nom ~ 0.588)
    fb = logspace(log10(cfg.m7_binding_band(1)), log10(cfg.m7_binding_band(2)), 2000).';
    b_nom = osi_ncf_looponly(squeeze(freqresp(ctx.C*ctx.H0d, 2*pi*fb)));

    buf = cfg.m7_buf;  Nbuf = round(buf*fs);
    nfft = 2^nextpow2(round(2*fs));  seg = min(2048, Nbuf);  ov = round(seg/2);
    win = hann(seg);
    tend = (cfg.t_start + buf) : cfg.m7_hop : cfg.t_end;
    nW = numel(tend);  tc = tend.' - buf/2;

    Spk = nan(nW,1);
    for k = 1:nW
        i1 = find(ctx.t >= tend(k) - buf, 1);
        y  = detrend(ctx.dm(i1:i1+Nbuf-1));          % OUTPUT ONLY
        [pp, ff] = pwelch(y, win, ov, nfft, fs);
        peak = max(   pp(ff>=cfg.m7_peak_band(1)  & ff<=cfg.m7_peak_band(2)));
        flr  = median(pp(ff>=cfg.m7_floor_band(1) & ff<=cfg.m7_floor_band(2)));
        Spk(k) = sqrt(peak / max(flr, eps));         % ~ ||S||_inf
    end
    Spk = medfilt1(Spk, cfg.m7_medfilt);

    Sref = median(Spk(tc < cfg.t_start + cfg.m7_calib_s), 'omitnan');   % calibrate at start
    ncf  = min(b_nom * Sref ./ Spk, cfg.m7_cap);                        % b(t) = b_nom*Sref/Spk

    below = ncf < cfg.ncf_threshold;
    hit = find(movsum(below, [cfg.m7_detect_consec-1, 0]) >= cfg.m7_detect_consec, 1);
    detect_time = NaN; if ~isempty(hit), detect_time = tc(hit); end

    out = struct('tc', tc, 'ncf', ncf, 'nominal', b_nom, 'threshold', cfg.ncf_threshold, ...
        'detect_time', detect_time, 'Spk', Spk, 'Sref', Sref);
end
```

Notes on the dependencies:
- `ctx.C`, `ctx.H0d`, `ctx.dm`, `ctx.t`, `ctx.fs` come from `osi_load_context`.
- `osi_ncf_looponly(L)` = `min |1+L| / sqrt(2(1+|L|²))` (the repo’s loop-only NCF).
- Only `dm` is used — no `dc`, no `dd`, no model FRF estimate.

## 9. How to build it from scratch (checklist)

1. **Add the config block** (`cfg.m7_*`) to `osi_config.m` (Section 7 table).
2. **Write `osi_method7_resonance.m`** in `common/` (Section 8).
3. **Write the runner** `trial_7/run_trial_7.m`:
   ```matlab
   this_dir = fileparts(mfilename('fullpath'));
   addpath(fullfile(this_dir,'..','common'));
   cfg = osi_config();  ctx = osi_load_context(cfg);
   t0 = tic;  m7 = osi_method7_resonance(ctx, cfg);  runtime = toc(t0);
   result = struct('tag','trial_7', ...
       'method', sprintf('Method 7 - Output resonance / sensitivity-peak (%.0fs buf)', cfg.m7_buf), ...
       'tc', m7.tc, 'ncf', m7.ncf, 'nominal', m7.nominal, 'threshold', m7.threshold, ...
       'runtime', runtime, 'detect_time', m7.detect_time);
   osi_report_series(result, cfg);     % save .mat + .png + print
   ```
4. **Pick the bands** from a baseline FFT of `dm` in a healthy window: the peak
   band around the resonance, the floor band on clean noise away from line tones.
5. **Run on each data file** (set `setappdata(0,'osi_data_file',path)` first), then
   `show_all_indicators.m` / `run_all_and_compare.m` pick it up automatically.

## 10. Validation (both records, default config)

| Record | `b_nom` | `Spk_ref` | stable-phase `b` | detection | end behaviour |
|---|---|---|---|---|---|
| **1731** broadband | 0.588 | ~2.3 | ~0.44–0.56 (no false alarm) | **717 s** (FRF methods: 728–734 s) | collapses with the loop |
| **1749** low-freq | 0.588 | ~1.5 | ~0.48–0.60 (no false alarm) | **219 s**; `b → 0.07` at 256–270 s | recovers (controller catches it) |

On 1749 the FRF methods (M1–M6) give **no detection** (they are blind at 35 Hz);
Method 7 is the only one that fires. On 1731 it cross-validates them.

## 11. Extensions

- **Closed-loop pole damping (alternative read-out).** Band-pass `dm` to 25–45 Hz,
  fit a low-order AR model, take the dominant pole’s damping `ζ`; `ζ → 0` at
  instability. More physical but noisier than the sensitivity peak here.
- **Adaptive bands.** Track the resonance centre (it shifts up near divergence) and
  move the peak band with it.
- **Saturation guard.** Pair with a raw band-energy trip for the fully-diverged
  regime where the peak/floor ratio saturates.
- **Method 8 — dither to recover the FRF.** Inject a small broadband probing
  signal (a dither at ~35 Hz) into the reference. That breaks the `dc = −C·dm`
  lock, the coherence drops below 1, and the FRF methods recover the *true*
  growing `H_dc_dm(35)` directly — turning the unobservable case into the
  observable (1731-like) case. This is the principled way to make `H_dc_dm`
  itself usable under a low-frequency operating input.