Sea regularity brainstorming

[mgrouch]

Great idea. You already maintain the machinery (I/Q envelope, P_{\text{disp}}=|\hat\eta|^2, and spectral moments). Here are drop-in metrics you can estimate online, with compact formulas and how they relate to “regularity”.

⸻

Moment-family (add M_3, M_4)

Let M_k=\langle P_{\text{disp}};\omega^k\rangle (your EMAs). You already have M_0,M_1,M_2.
• Add: M_3=\langle P\omega^3\rangle,; M_4=\langle P\omega^4\rangle.

From these:
• Mean angular frequency: \bar\omega = M_1/M_0 (you already use this).
• Central moments (about \bar\omega):
\mu_2 = \frac{M_2}{M_0}-\bar\omega^2,\qquad
\mu_3 = \frac{M_3}{M_0}-3\bar\omega\frac{M_2}{M_0}+2\bar\omega^3,
\mu_4 = \frac{M_4}{M_0}-4\bar\omega\frac{M_3}{M_0}+6\bar\omega^2\frac{M_2}{M_0}-3\bar\omega^4.
• (Excess) kurtosis (your “M4 metric”):
\gamma_2 = \frac{\mu_4}{\mu_2^2}-3.
Interpretation: narrow, peaky spectra \Rightarrow \gamma_2 \uparrow. Broadband \Rightarrow \gamma_2 \downarrow (toward 0 for quasi-Gaussian).
• Skewness:
\gamma_1 = \frac{\mu_3}{\mu_2^{3/2}}.
Interpretation: asymmetric spectral mass around \bar\omega (e.g., wind-sea tail vs swell).
• Classical oceanographic bandwidth (Cartwright–Longuet-Higgins):
\varepsilon = \sqrt{1-\frac{M_2^2}{M_0,M_4}} \quad \in [0,1].
Interpretation: \varepsilon!\downarrow \Rightarrow more narrowband (regular); \varepsilon!\uparrow broadband.
• Period summaries (very handy readouts):
T_{m01}=\frac{2\pi,M_0}{M_1},\qquad
T_{m02}=2\pi\sqrt{\frac{M_0}{M_2}}.
• Mean zero-upcrossing period: T_z = T_{m02}.
• Upcrossing rate: \nu_0=\dfrac{1}{T_z}=\dfrac{1}{2\pi}\sqrt{M_2/M_0}.

⸻

Entropy-family (broadband vs regular)

Let a power-weighted frequency PDF be
p(\omega) \approx \frac{P_{\text{disp}}(\omega)}{\int P_{\text{disp}}(\omega’),d\omega’}.
Online, approximate by binning \omega (e.g., 16–32 bins) and EMA-accumulating E_i=\langle P\cdot \mathbf{1}{\omega\in \text{bin }i}\rangle, then p_i=E_i/\sum_j E_j.
• Shannon spectral entropy (discrete):
H{\text{spec}} = -\sum_i p_i \log p_i,\qquad
H_{\text{norm}}=\frac{H_{\text{spec}}}{\log N_{\text{bins}}}\in[0,1].
Interpretation: regular/narrow \Rightarrow H_{\text{norm}}!\downarrow; broadband \Rightarrow H_{\text{norm}}!\uparrow.
• Rényi entropy (order \alpha, e.g., 2):
H_\alpha = \frac{1}{1-\alpha}\log!\sum_i p_i^\alpha
More sensitive to peaks for \alpha>1.
• Spectral flatness (geometric/arith. mean):
\mathrm{SF} = \frac{\exp\big(\sum_i p_i \log s_i\big)}{\sum_i p_i s_i},\quad s_i:=E_i/(\Delta\omega_i)
(Or compute on envelope power samples with a log-EMA and linear-EMA pair.)
Interpretation: \mathrm{SF}\to 0 for peaky (regular), \mathrm{SF}\to 1 for white-like.

⸻

Phase/circular statistics add-ons

You already track R_{\text{phase}} = |\langle \mathbf{u}\rangle| with \mathbf{u} = z/|z|.
• Von Mises concentration \kappa from R (use standard R!\to!\kappa approximation); then circular (differential) entropy of the von Mises model:
h_{\text{circ}} \approx \log!\big(2\pi I_0(\kappa)\big) - \kappa,\frac{I_1(\kappa)}{I_0(\kappa)}.
Interpretation: lower h_{\text{circ}} \Rightarrow more phase-locked (regular).
(If you’d rather avoid Bessel: coarse-bin the phase of z into p_k and compute Shannon entropy like above.)
• Phase slip rate: EMA of |\dot\phi - \omega_{\text{lp}}| normalized by \omega_{\text{lp}}. Interpretation: demodulation consistency; lower is more regular.
• Rice K-factor on the complex envelope z:
K = \frac{|\langle z\rangle|^2}{\langle |z-\langle z\rangle|^2\rangle}.
Interpretation: large K = strong “line-of-sight” (single tone) + weak scatter; small K = diffuse/broadband.

⸻

Amplitude/envelope statistics

Using |\hat\eta| (disp envelope magnitude):
• Crest factor:
\mathrm{CF} = \frac{\max_t |\hat\eta|}{\sqrt{\langle \hat\eta^2\rangle}}\quad
\text{(use running max over, e.g., (30–60) s).}
Higher CF often indicates intermittent broadening or harmonics.
• Envelope coefficient of variation:
\mathrm{CV} = \frac{\mathrm{std}(|\hat\eta|)}{\mathrm{mean}(|\hat\eta|)}.
Regular seas trend to lower CV (after spin-up).
• Rayleighness test (for narrowband Gaussian seas): compare empirical CDF of |\hat\eta| to Rayleigh with parameter \sigma^2=\tfrac{1}{2}\langle \hat\eta^2\rangle; report a KS distance D_R. Small D_R supports classical narrowband assumption.

⸻

Quick oceanographic “derived” set (from M_0,M_1,M_2,M_4)
• Significant height (you already have): H_s\approx 4\sqrt{m_0} for displacement; your online M_0 is the time-domain analog—keep your blend.
• Peak period proxies:
• T_{m01}=2\pi M_0/M_1
• T_{m02}=2\pi\sqrt{M_0/M_2}=T_z (mean zero-upcrossing period)
• Bandwidth \varepsilon as above, and optionally the Benassai parameter B = \sqrt{1-\dfrac{M_1^2}{M_0 M_2}} (another width measure, less sensitive to tails).

⸻

Implementation sketch (what to add to your class)
• Accumulators: M3, M4 (float), optional bins_K (small array) for entropy, plus EMAs for meanZ and varZ to get Rice K.
• In updateSpectralMoments():

// After computing P_disp and omega_norm
M3 = (1-a_m)M3 + a_m * (P_disp * omega_normomega_normomega_norm);
M4 = (1-a_m)M4 + a_m * (P_disp * omega_normomega_normomega_norm*omega_norm);

•	Maintain meanZ_r/i and varZ over the displacement envelope (or over z directly) to compute K.
•	For entropy: choose, say, 24 ω-bins between [\omega_{\min},\omega_{\max}]; EMA add P_{\text{disp}} into the active bin; readout H_\text{norm} from normalized bin sums.

⸻

How they help you
• Diagnostics: \gamma_2 (excess kurtosis) and \varepsilon separate “really narrow” (clean swell) from “merely not too wide”; entropy and flatness catch multi-peak/shouldered spectra your RBW can miss.
• Fusion: you can build a richer R by blending R_{\text{phase}} with two orthogonal “spread” signals: one moment-based (\varepsilon or \gamma_2) and one information-theoretic (H_{\text{norm}} or SF).
• QA: T_{m01}, T_{m02} give period sanity checks vs your tracker; K exposes hidden AM/FM clutter.

If you want, I can patch your header with M3/M4 + ε + γ2 + H_norm + K fields and exact lines to compute each in your existing methods.

[mgrouch]

My frequency estimates are noisy and signal is noisy too with some constant bias. How accurate are these metrics under these conditions

[mgrouch]

Short answer: most of your metrics still work, but they bias “too messy” when (i) the frequency estimate jitters and (ii) the accel signal carries DC/low-freq bias. Below is a quick “how accurate / how it fails / how to harden” guide for each metric, with first-order error scalings so you can judge impact.

⸻

What frequency/accel noise does to your pipeline

Let \omega=\bar\omega+\delta with \mathbb{E}[\delta]=0, \mathrm{Var}(\delta)=\sigma_\omega^2. You form displacement envelope via \hat\eta \sim z/\omega^2, hence powers \propto 1/\omega^4.
• Bias of inverse powers (Jensen):
\mathbb{E}!\left[\frac{1}{\omega^2}\right] \approx \frac{1}{\bar\omega^2}\Big(1+ \frac{3\sigma_\omega^2}{\bar\omega^2}\Big),\quad
\mathbb{E}!\left[\frac{1}{\omega^4}\right] \approx \frac{1}{\bar\omega^4}\Big(1+ \frac{10\sigma_\omega^2}{\bar\omega^2}\Big).
So M0 (power) tends to be biased high, but… frequency moments and RBW can be biased too (see below).
• DC/low-freq accel bias b: before demod it’s constant; after demod it becomes b,e^{-j\phi(t)} (a pure tone at \omega). Your envelope EMA largely rejects it if its bandwidth \ll \omega, but if \alpha_\text{env} is high (short \tau_\text{env}), some leakage remains.

⸻

Metric-by-metric accuracy

R_{\text{phase}} (mean resultant length)
• Sensitivity: Only to phase jitter of the demod carrier. If phase error has small circular variance \sigma_\theta^2, then
R_{\text{phase}} \approx 1 - \tfrac{1}{2}\sigma_\theta^2.
• Robust to: amplitude noise, DC bias (it uses unit vectors).
• Degrades with: noisy \omega → noisy \phi → larger \sigma_\theta^2.
• Harden: keep your slew-limited \omega; increase \tau_\text{coh} slightly; apply the tiny-magnitude fade (prevents sticky high R when |z| is tiny).

Bottom line: high accuracy if the tracker is decent; it’s your most robust regularity channel.

⸻

R_{\text{spec}} via RBW (two-pole \mu_2)
• Sensitivity: to \omega jitter directly; your subtraction of \mathrm{Var}(\omega) helps but isn’t perfect if noise statistics shift.
• Bias direction: jitter inflates RBW ⇒ R_{\text{spec}} biased low (appears less regular).
• Harden: what you’re already doing:
• subtract slow variance (\text{var_slow}) and clamp pre-sqrt,
• lengthen \tau_\text{mom} a bit,
• optionally compute a parallel RBW from phase increments (variance of \Delta\phi/dt); fuse the two by min/geomean.

Bottom line: moderately accurate; expect conservative (low) bias when tracker is noisy.

⸻

Moment family: M_0\ldots M_4, \mu_2,\mu_3,\mu_4, skew/kurtosis, bandwidth \varepsilon
• Sensitivity: strong to \omega jitter because P_{\text{disp}}\propto 1/\omega^4.
• Bias: M_k inflated by \sim (1+c_k,\sigma_\omega^2/\bar\omega^2) with c_0\approx10 for the 1/\omega^4 part; ratios like \bar\omega=M_1/M_0 partly cancel but not fully (cross-terms remain).
• Kurtosis/excess \gamma_2=\mu_4/\mu_2^2-3: noisy \omega tends to push \gamma_2 up (looks peakier) because \mu_4 grows faster than \mu_2^2.
• Cartwright–Longuet-Higgins bandwidth \varepsilon=\sqrt{1-M_2^2/(M_0 M_4)}: frequency jitter pushes \varepsilon up (looks more broadband).
• Harden:
• apply your 1st-order 1/\omega^2 bias correction (and extend to 1/\omega^4 by reusing the same factor squared if you want),
• increase \tau_\text{omega} slightly so \omega_\text{norm} is smoother than envelope dynamics,
• compute parallel moments on raw accel envelope z (no 1/\omega^2) just for shape (skew/kurt/entropy), then map conclusions back to displacement.

Bottom line: usable trends but numerically touchy; treat as supporting evidence unless \omega is clean.

⸻

Entropy / flatness on spectral bins
• Sensitivity: very sensitive to binning vs. \omega jitter: jitter smears p_i across bins ⇒ entropy is biased high, flatness increases.
• Harden:
• choose bin width \Delta\omega \gtrsim 2\sigma_\omega (estimate \sigma_\omega from your slow variance),
• use EMA bins with the same \tau_\text{mom},
• apply a small Miller–Madow correction for discrete entropy bias if sample counts are low (optional),
• or compute entropy on phase-increment PDF instead of \omega directly (often cleaner).

Bottom line: reliable directionally if you right-size bins; absolute numbers will be conservative (higher H).

⸻

Rice K on complex envelope
• If computed on z (accel envelope):
• Insensitive to \omega noise (no 1/\omega^2 factor),
• Sensitive to additive accel noise (raises denominator) ⇒ K biased low.
• If computed on displacement envelope z/\omega^2: inherits all 1/\omega^2 bias problems.
• Harden: compute K on z, not on displacement, and optionally whiten z with a short AR(1) or a fixed pre-LP.

Bottom line: good robustness if you keep it in the accel domain.

⸻

H_s estimate
• Sensitivity: to both \omega jitter (through 1/\omega^4) and additive noise (through |\hat\eta|^2).
• Bias direction: with your current mapping, regular seas tend to be underestimated if \omega noise depresses R and if you under-correct 1/\omega^2.
• Harden: keep the 1/\omega^2 bias correction, use R-aware blending (you already do), and consider a small noise floor subtraction on P_{\text{disp}} estimated from calm intervals.

Bottom line: accurate to ≈±10–20% in practice after corrections; without them it can be off more.

⸻

Practical “accuracy tier” summary (with noisy \omega and biased accel)
• High-robustness: R_{\text{phase}}, Rice K on z, zero-upcrossing period T_{m02} (if computed from phase increments rather than \omega).
• Medium: RBW-based R_{\text{spec}} (with variance subtraction), entropy/flatness (with bin size ≥ jitter).
• Low (treat as diagnostic only unless \omega is clean): kurtosis/skewness from M_3,M_4, \varepsilon bandwidth, displacement-domain K.

⸻

Concrete hardening you can ship (cheap and effective)
1. DC/low-freq bias remover on accel (before demod):
EMA with \tau_\text{bias}=60\text{–}120 s:

bias = (1-a_b)bias + a_baccel_z; accel_z -= bias;

2.	Bias correction for 1/\omega^2 (you already have the 1st-order form):

Use it for powers, too (i.e., square the correction when forming P_{\text{disp}}).
3. Phase-increment track for a “clean” RBW:
In parallel, set \tilde\omega = \Delta\phi/dt (from your slewed \phi), keep its mean/var EMAs, and compute an RBW’ independent of the external \omega_\text{inst}. Fuse with your existing RBW using geometric mean or minimum.
4. Entropy bins matched to jitter:
Set \Delta\omega = 2\sqrt{\text{var_slow}} (cap between, say, 0.02–0.15 rad/s), 16–24 bins.
5. Report CIs / stability flags:
Expose a quality bit based on \sigma_\omega^2 / \bar\omega^2 and envelope SNR; when “poor,” down-weight entropy/kurtosis in your fused R.

⸻

If you want, I can annotate your header with the exact extra lines (bias remover, squared correction for P_{\text{disp}}, and a tiny phase-RBW block) so you can toggle them via #define and see the deltas quickly.

[mgrouch]

And for wave height envelope

[mgrouch]

Got it—wave height envelope specifically. Here’s how to make your H_s (and the instantaneous height envelope) accurate when ω is noisy and the acceleration has DC/low-freq bias.

⸻

What to fix (in priority order)

1. Remove accel bias before anything

A slow EMA on a_z avoids baking DC into the demod:

// member vars
float a_bias = 0.0f;
constexpr static float TAU_BIAS = 90.0f; // 1–2 min
float alpha_bias = 0.0f;

// in updateAlpha(dt)
alpha_bias = 1.0f - std::exp(-dt_s / TAU_BIAS);

// in update(...)
a_bias = (1.0f - alpha_bias) * a_bias + alpha_bias * accel_z;
accel_z -= a_bias;

2. Correct the 1/ω² bias (you already have 1st-order—apply it to power too)

Because P_{\text{disp}} \propto 1/\omega^4, square the correction used for 1/\omega^2:

// after computing omega_norm, var_slow
float inv_w2_den = std::max(omega_norm * omega_norm, EPSILON);
float corr = 1.0f + std::clamp(var_slow / inv_w2_den, 0.0f, 0.5f);
float inv_w2 = corr / inv_w2_den;

// envelope → displacement
float disp_r = z_real * inv_w2;
float disp_i = z_imag * inv_w2;
float P_disp = disp_rdisp_r + disp_idisp_i;

// apply squared correction to power implicitly via inv_w2^2
// (done already since disp uses inv_w2)

This alone typically lifts underestimates in regular seas by 5–15%.

3. Subtract a noise floor from P_{\text{disp}} before moments

Track a conservative floor and clamp:

// member vars
float P_floor_fast = 0.0f, P_floor_slow = 0.0f;

// in updateSpectralMoments(), after P_disp
P_floor_fast = (1.0f - 0.10f) * P_floor_fast + 0.10f * P_disp;
P_floor_slow = (1.0f - 0.01f) * P_floor_slow + 0.01f * P_disp;
float Pn = std::min(P_floor_fast, P_floor_slow); // min-EMA floor tracker
float P_corr = std::max(0.0f, P_disp - 0.5f * Pn); // subtract half the floor (stable)

// use P_corr (not P_disp) for M0..M2 (and M3/M4 if you add them)

This shaves off sensor/electrical noise and residual demod leakage.

4. Make ω-blend more phase-aware

You already blend omega_lp with omega_disp_lp. Use convex mapping:

float lambda = LAMBDA_MIN + (LAMBDA_MAX - LAMBDA_MIN) * (R_phase * R_phase);
lambda = std::clamp(lambda, LAMBDA_MIN, LAMBDA_MAX);
float omega_norm = has_omega_disp_lp
? lambda * omega_lp + (1.0f - lambda) * omega_disp_lp
: omega_lp;

This reduces 1/ω² wobble when phase coherence is mediocre.

5. Use an R-aware H_s blend and a tiny “regular-sea lift”

Keep your H_s \approx C_H \cdot 2\sqrt{M_0}\cdot f(R). Improve two parts:
• R-aware smoother: you already adapted R_{\text{out}} smoothing—good.
• Lift for R>0.7: a 1–3% gain counters conservative bias once seas are clearly regular:

static float heightFactorFromR(float R) {
float x = std::clamp(R, 0.0f, 1.0f);
float s = xx(3.0f - 2.0f*x); // smoothstep
const float C_reg = 1.0f, C_irr = std::sqrt(2.0f);
float base = C_irr - (C_irr - C_reg) * s;
float tilt = 0.10f * std::clamp((x - 0.7f) / 0.3f, 0.0f, 1.0f); // up to -10%
return base * (1.0f - 0.12f * tilt); // ~1–3% lift for very regular
}

⸻

Instantaneous height envelope (crest-to-trough proxy)

Sometimes you want a sample-level height envelope H_{\text{env}}(t), not just H_s. Do it on the displacement envelope magnitude r=|\hat\eta|:
1. Compute a robust magnitude with log-smoothing (reduces spikes):

// member vars
float log_r_lp = 0.0f; bool has_log=false;

// per update, after disp_r/disp_i
float r = std::sqrt(disp_rdisp_r + disp_idisp_i);
float log_r = std::log(std::max(r, 1e-9f));
if (!has_log) { log_r_lp = log_r; has_log = true; }
else { log_r_lp = (1.0f - alpha_env) * log_r_lp + alpha_env * log_r; }

float r_smooth = std::exp(log_r_lp);

2.	Map to height envelope with the same f(R) you use for H_s:

float H_env_inst = 2.0f * r_smooth * heightFactorFromR(R_out) * H_CAL;

3.	Stabilize with a “crest-tracker” pair: a slow EMA for the median and a faster max-tracker for peaks; combine to avoid wind-gust spikes:

// member vars
float H_med = 0.0f, H_peak = 0.0f;
float a_med = 1.0f - std::exp(-last_dt / 20.0f); // ~20 s
float a_pk = 1.0f - std::exp(-last_dt / 5.0f); // ~5 s

H_med = (1.0f - a_med) * H_med + a_med * H_env_inst;
H_peak = std::max((1.0f - a_pk) * H_peak, H_env_inst); // leaky peak
float H_env = 0.65f * H_peak + 0.35f * H_med; // robust envelope readout

This gives you a stable, responsive height envelope that tracks crests but ignores transient noise.

⸻

Sanity/QA checks you can expose
• Quality bit: if \sigma_\omega^2/\bar\omega^2 > 0.02 (e.g., >14% std/mean), flag Low confidence; de-emphasize H_{\text{env}} changes.
• Floor dominance: if P_n/P_{\text{corr}} > 0.3, you’re noise-limited—report H_s with caution.

⸻

Expected accuracy after these fixes
• Regular seas (good coherence, modest ω jitter): H_s within ±10–15%, envelope tracks crests smoothly without “jerks”.
• Messy seas (low R, high jitter): conservative bias remains by design (safer to under-report a bit). You can still keep ±20–25% if P_n isn’t dominant.

If you want, I can splice these blocks directly into your header (minimal diffs around updateSpectralMoments(), computeRegularityOutput(), and a tiny accessor getWaveHeightEnvelopeInstant()).