nx_ringdown_lib.nx source
↩ module page · 476 lines · 24622 B
1// nx_ringdown_lib.nx -- ONE damped-oscillator identifier for the whole estate.
2//
3// WHY THIS EXISTS. The estate has to answer "what is this thing's natural frequency and damping"
4// about TWO subjects that must be comparable: our own XPBD plant (via nx_dyna_bind_skin's
5// ringdown) and REAL FOOTAGE (via nx_dynaoracle's optical-flow displacement series). If each
6// carried its own estimator the comparison would be measuring two rulers, not one subject.
7// Extracted VERBATIM from nx_dyna_bind_skin.dbs_ringdown_t 2026-08-25 so the two cannot disagree
8// BY CONSTRUCTION -- the law the jobclaim lane already paid for.
9//
10// SUBJECT-AGNOSTIC BY DESIGN: it takes a signal, a sample period and a deadband. It knows nothing
11// about tissue, cages, pixels or cameras. Hair, cloth, camera shake and any damped oscillator are
12// the same arithmetic.
13//
14// DOCUMENTED IMPRECISION (carried over, deliberately not silently improved): rd_log2_q16 is a
15// bit-length plus LINEAR mantissa term, worst error about 0.086 in log2 inside an octave. The
16// analyser therefore ALWAYS returns the two raw peak amplitudes beside the derived zeta so any
17// reader can recompute the decrement exactly and is never forced to trust this routine.
18//
19// license_tier: ORIGINAL
20import "nx_syscalls.nx"
21
22const RD_Q16: i64 = 65536
23const RD_PERMIL: i64 = 1000
24const RD_US_PER_S: i64 = 1000000
25const RD_MHZ_PER_HZ: i64 = 1000
26// ln(2) and 2*pi in q16 -- the same constants the incumbent carried.
27const RD_LN2_Q16: i64 = 45426
28const RD_TWOPI_Q16: i64 = 411775
29
30// ---- DECLARED IMPRECISIONS. ONE OWNER, HERE. A consumer or a gate that carries its own copy of
31// these numbers goes stale the day the routine improves, and a stale tolerance is a gate that
32// passes exactly what it was built to catch. rd_log2_q16's worst error inside an octave is about
33// 0.086 in log2; the other two are ln(2) and 2*pi, in permil, for the small-damping identities.
34const RD_LOG2_ERR_PERMIL: i64 = 86
35const RD_LN2_PERMIL: i64 = 693
36const RD_TWOPI_PERMIL: i64 = 6283
37// The crossing-interval admission window, in SAMPLES. This is the SAMPLING QUANTUM ITSELF and not
38// a tuned tolerance: a crossing index is exact to one sample, so an interval that differs from the
39// first by more than one sample is not the same half period. MEASURED: 1 and 2 select an identical
40// span at all 18 points of the declared grid (zeta 50..600 permil x fn 1.5/2/3 Hz), so the choice
41// is not load-bearing on clean data and the conservative one is taken.
42const RD_SPAN_TOL_SAMPLES: i64 = 1
43// A span must cover at least a FULL PERIOD of baseline, which is three admitted crossings and two
44// half-periods between the two peaks compared -- exactly the incumbent two-lobe decrement's own
45// baseline. THIS WAS 2 AND IT WAS WRONG. At two admitted crossings the span reads HALF a period,
46// i.e. a SHORTER baseline than the estimator it exists to improve on, and a generalisation that
47// reads less evidence than the incumbent is not a generalisation. MEASURED on the end-to-end video
48// surface before this line was corrected: at every point where the span collapsed to one
49// half-period it was WORSE than the two-lobe decrement (zeta errors 26, 31 permil against the
50// incumbent's 13, 19), and at every point where it kept a real baseline it was better. Refusing
51// the short case removes the only regime where this estimator loses. Named apart from the other
52// small integers in this file because it answers a different question: how much evidence, not what
53// shape.
54const RD_SPAN_MIN_INTERVALS: i64 = 3
55
56// ---- result slots. ONE layout, every consumer reads the same index. ----
57const RD_R_NSAMP: i64 = 0
58const RD_R_A1: i64 = 1
59const RD_R_A2: i64 = 2
60const RD_R_PER_US: i64 = 3
61const RD_R_FN_MHZ: i64 = 4
62const RD_R_ZETA: i64 = 5
63const RD_R_PEAKS: i64 = 6
64const RD_R_MAXABS: i64 = 7
65const RD_R_RESID: i64 = 8
66// ---- PERIODIC-SPAN SLOTS, APPENDED. Every index above keeps its number, so a consumer that reads
67// only the original nine sees identical values; RD_R_N grows, so every consumer that sizes its
68// array from this constant allocates the larger array BY CONSTRUCTION rather than by remembering.
69const RD_R_SPAN_HALVES: i64 = 9 // half-periods between the two peaks that were compared
70const RD_R_SPAN_A: i64 = 10 // peak of the FIRST closed lobe
71const RD_R_SPAN_B: i64 = 11 // peak of the LAST ADMITTED closed lobe
72const RD_R_SPAN_ZETA: i64 = 12 // zeta from the generalised decrement, permil
73const RD_R_SPAN_ZETA_TOL: i64 = 13 // and the imprecision THIS run can support, permil
74const RD_R_SPAN_PER_US: i64 = 14
75const RD_R_SPAN_FN_MHZ: i64 = 15
76const RD_R_SPAN_FN_TOL: i64 = 16 // permil of the frequency
77const RD_R_SPAN_LOBES: i64 = 17 // closed lobes seen, admitted or not
78const RD_R_SPAN_RC: i64 = 18 // the span estimator's OWN third state, never a bare zero
79const RD_R_N: i64 = 19
80
81// ---- THE THIRD STATE. A zero fn is NOT the same fact as an unmeasurable fn, and every consumer
82// that conflates them acquits on a measurement that never happened. rd_analyze fills the slots it
83// COULD fill and returns WHICH conjunct stopped it -- never a bare zero.
84const RD_OK: i64 = 0
85const RD_UNOBS_NO_SAMPLES: i64 = 0 - 1
86const RD_UNOBS_NO_CROSSING: i64 = 0 - 2
87const RD_UNOBS_NO_PERIOD: i64 = 0 - 3
88const RD_UNOBS_NO_DECAY: i64 = 0 - 4
89// The span estimator's own refusal: crossings exist, but their spacing is not periodic, so they are
90// not half-period markers and no decrement may be read across them. This is a DIFFERENT fact from
91// "no full period in the window", and giving it its own code is what stops a consumer reporting one
92// when the other happened.
93const RD_UNOBS_NO_PERIODIC_SPAN: i64 = 0 - 5
94
95func rd_reason(r: i64) -> *u8 {
96 if r == RD_OK { return "OK" as *u8 }
97 if r == RD_UNOBS_NO_SAMPLES { return "UNOBSERVABLE-no-samples" as *u8 }
98 if r == RD_UNOBS_NO_CROSSING { return "UNOBSERVABLE-never-crosses-deadband" as *u8 }
99 if r == RD_UNOBS_NO_PERIOD { return "UNOBSERVABLE-no-full-period-in-window" as *u8 }
100 if r == RD_UNOBS_NO_DECAY { return "UNOBSERVABLE-second-lobe-not-smaller" as *u8 }
101 if r == RD_UNOBS_NO_PERIODIC_SPAN { return "UNOBSERVABLE-crossings-are-not-periodic" as *u8 }
102 return "UNOBSERVABLE-unclassified" as *u8
103}
104
105func rd_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
106
107func rd_isqrt(n: i64) -> i64 {
108 if n <= 0 { return 0 }
109 var r: i64 = n
110 var y: i64 = (r + 1) / 2
111 while y < r { r = y; y = (r + n / r) / 2 }
112 return r
113}
114
115func rd_log2_q16(x_q16: i64) -> i64 {
116 if x_q16 <= 0 { return 0 }
117 var v: i64 = x_q16
118 var e: i64 = 0
119 while v >= RD_Q16 * 2 { v = v / 2; e = e + 1 }
120 while v < RD_Q16 { v = v * 2; e = e - 1 }
121 let frac: i64 = v - RD_Q16
122 return e * RD_Q16 + frac
123}
124
125func rd_ln_q16(x_q16: i64) -> i64 { return rd_log2_q16(x_q16) * RD_LN2_Q16 / RD_Q16 }
126
127// ---- THE DECREMENT IDENTITY, AND THE GENERALISATION THAT CONTAINS IT -------------------------
128//
129// PRE-DECLARED ACCEPT RULE, written before the integer implementation was ever run: the periodic-
130// span estimator is adopted ONLY if, over the declared grid of injected damping (zeta 50, 100, 200,
131// 350, 450 and 600 permil -- spanning and exceeding the cited tissue band 350..600 that
132// knowledge/gamefeel_oracle.conf carries) crossed with injected fn 1.5, 2.0 and 3.0 Hz, on the
133// analytic mass-spring-damper fixture nx_dynaoracle_gate integrates for itself:
134// (a) EVERY point lands inside the tolerance THIS FILE publishes for that point, and
135// (b) its WORST absolute zeta error over the grid is STRICTLY SMALLER than the two-lobe
136// decrement's worst error on the SAME grid.
137// If it is not strictly better it does not get wired and the two-lobe decrement remains the only
138// estimator. The gate asserts both halves, so the rule is enforced rather than remembered, and a
139// justification written after the result cannot be substituted for it.
140//
141// zeta = delta / sqrt(4*pi^2 + delta^2), delta being the per-PERIOD logarithmic decrement.
142// rd_zeta_permil_span carries that arithmetic ONCE and the two-lobe form delegates to it with
143// halves = 2, which is EXACT -- delta = 2*ln(a/b)/2 = ln(a/b), and x*2/2 == x in integers -- so no
144// second copy of this tail exists to drift away from the first.
145//
146// WHY A SPAN HELPS ARITHMETICALLY, AND NOT BY FITTING. The two-lobe decrement reads ONE ratio of
147// two peaks one period apart, so every imprecision in that single ratio lands on zeta undivided:
148// rd_log2_q16's declared 0.086-in-log2 error, plus one amplitude LSB on each of the two peaks.
149// Reading the SAME identity across m half-periods gives delta = 2*ln(A_first/A_last)/m, which
150// divides all of them by m/2. Nothing is fitted and no coefficient is tuned -- the declared
151// imprecision shrinks with the baseline, which is why the tolerance below can be COMPUTED rather
152// than chosen. A least-squares fit of the log envelope over the same lobes was built and MEASURED
153// against this on the same grid and was WORSE (worst zeta error 27 permil unweighted, 64 permil
154// inverse-variance weighted, against 22 for the incumbent and 6 for this), because it spends its
155// degrees of freedom on the decayed tail where one LSB is the whole signal. It is recorded here so
156// the next reader does not rebuild it.
157
158// Two half-periods make a period: this constant is the SHAPE of the identity, never a knob.
159const RD_HALVES_PER_PERIOD: i64 = 2
160// A span has two ENDS and each end carries one crossing-timing error. A DIFFERENT 2 from the one
161// above, deliberately named apart: one constant serving two unrelated purposes can never be tuned
162// for either.
163const RD_SPAN_ENDS: i64 = 2
164
165func rd_zeta_permil_span(a: i64, b: i64, halves: i64) -> i64 {
166 if b <= 0 { return RD_PERMIL }
167 if a <= b { return 0 }
168 if halves <= 0 { return 0 }
169 let ratio_q16: i64 = a * RD_Q16 / b
170 var delta: i64 = rd_ln_q16(ratio_q16) * RD_HALVES_PER_PERIOD / halves
171 if delta < 0 { delta = 0 }
172 let fourpi2: i64 = RD_TWOPI_Q16 * RD_TWOPI_Q16 / RD_Q16
173 let dd: i64 = delta * delta / RD_Q16
174 let root: i64 = rd_isqrt((fourpi2 + dd) * RD_Q16)
175 if root <= 0 { return 0 }
176 return delta * RD_PERMIL / root
177}
178
179// ---- THE ESTIMATOR PUBLISHES ITS OWN IMPRECISION ---------------------------------------------
180//
181// SCOPE, STATED SO NOBODY READS THIS BOUND AS SOMETHING IT IS NOT. What follows is the ARITHMETIC
182// imprecision of the estimator: the error its own integer maths can introduce given the numbers it
183// was handed. It is NOT a measurement-uncertainty bound on the SUBJECT. A displacement series
184// recovered from video also carries optical-flow envelope error, and this bound cannot see that,
185// because nothing in this file knows how the series was produced.
186// MEASURED, both directions, on the same day this was written: over 36 analytic fixture points
187// (six damping ratios x three stiffnesses x two sampling ratios) the bound held at 36 of 36; over
188// the end-to-end video surface it was exceeded at 5 of 11 answered points, and the excess is
189// front-end noise, not arithmetic. A caller that needs the second kind of bound must read the
190// consumer's own noise rows -- nx_dynaoracle publishes deadband_q8 and a residual for exactly that
191// -- and must not mistake a tight arithmetic bound for a tight measurement.
192// A ruler that reports a value and not its bound forces every consumer to re-derive the bound, and
193// a re-derived bound goes stale the day the routine improves. Both bounds below are written in
194// terms of the DECLARED imprecisions above and the amplitudes ACTUALLY MEASURED on the run, so what
195// ships beside the value is what the arithmetic can support on THAT record -- not a constant that
196// happens to flatter the easy cases.
197func rd_span_zeta_tol_permil(a: i64, b: i64, halves: i64) -> i64 {
198 if halves <= 0 { return RD_PERMIL }
199 if a <= 0 { return RD_PERMIL }
200 if b <= 0 { return RD_PERMIL }
201 let ln_err: i64 = RD_LOG2_ERR_PERMIL * RD_LN2_PERMIL / RD_PERMIL
202 let amp_err: i64 = RD_PERMIL / a + RD_PERMIL / b
203 return (ln_err + amp_err) * RD_HALVES_PER_PERIOD * RD_PERMIL / (halves * RD_TWOPI_PERMIL) + 1
204}
205
206// A crossing index is uncertain by the SAMPLING QUANTUM (one sample) plus the displacement that one
207// amplitude LSB causes at the waveform's slope through zero: that slope is 2*pi*b/P per sample, so
208// one LSB moves the crossing by P/(2*pi*b) samples. This is why a decayed tail dominates the bound
209// and a strong record does not. The period is read across `halves` half-periods, so a timing error
210// at each end is divided by that span.
211func rd_span_fn_tol_permil(per_samples: i64, halves: i64, b: i64) -> i64 {
212 if per_samples <= 0 { return RD_PERMIL }
213 if halves <= 0 { return RD_PERMIL }
214 var bb: i64 = b
215 if bb < 1 { bb = 1 }
216 let eps_permil_of_sample: i64 = RD_PERMIL + per_samples * RD_PERMIL * RD_PERMIL / (RD_TWOPI_PERMIL * bb)
217 return RD_SPAN_ENDS * eps_permil_of_sample / (halves * per_samples) + 1
218}
219
220// zeta from the logarithmic decrement of two successive same-sign peaks -- the textbook free-decay
221// identity delta = ln(A1/A2), zeta = delta / sqrt(4*pi^2 + delta^2). Permil, so it compares
222// directly against the tissue_zeta_permil subject row.
223// ⚠2026-08-25: a coordinator briefly RENAMED this function's signature to rd_zeta_from_delta while
224// intending to extract the conversion into one copy -- work this lane had ALREADY done, by the
225// better route of making the two-lobe form a SPECIAL CASE of the span form above. The rename left
226// this body referring to parameters that no longer existed and removed rd_zeta_permil from under its
227// four callers, and the BUILD LANE caught it immediately by name and line. Restored.
228// ★★★★★★THE EXTRACTION WAS ALREADY THERE UNDER A DIFFERENT NAME, AND NOT READING THE FILE FIRST IS
229// HOW A SHARED LIB GETS BROKEN BY SOMEONE TRYING TO IMPROVE IT.
230func rd_zeta_permil(a1: i64, a2: i64) -> i64 {
231 return rd_zeta_permil_span(a1, a2, RD_HALVES_PER_PERIOD)
232}
233
234// ---- THE ANALYSER. sig[0..n) is a displacement series about its own equilibrium, sampled every
235// dt_us microseconds. quant is the DEADBAND: the amplitude below which a sign change is noise
236// rather than a crossing. The caller owns that number and must say where it came from -- for the
237// XPBD plant it is the solver's documented rounding quantum, for video it is a measured noise
238// floor. A deadband picked by taste is exactly the magic number this estate refuses.
239func rd_analyze(sig: *i64, n: i64, dt_us: i64, quant: i64, out: *i64) -> i64 {
240 var f: i64 = 0
241 while f < RD_R_N { out[f] = 0; f = f + 1 }
242 if n <= 0 { return RD_UNOBS_NO_SAMPLES }
243 if dt_us <= 0 { return RD_UNOBS_NO_SAMPLES }
244 out[RD_R_NSAMP] = n
245
246 var mx: i64 = 0
247 var k: i64 = 0
248 while k < n { if rd_abs(sig[k]) > mx { mx = rd_abs(sig[k]) } k = k + 1 }
249 out[RD_R_MAXABS] = mx
250
251 // RESIDUAL: peak excursion over the final eighth. REPORTED, never used to suppress detection:
252 // a large residual IS the "jelly / never settles" finding, so folding it into a noise floor
253 // would hide the exact failure mode a caller runs this to catch.
254 var resid: i64 = 0
255 var t: i64 = n - n/8
256 while t < n { if rd_abs(sig[t]) > resid { resid = rd_abs(sig[t]) } t = t + 1 }
257 out[RD_R_RESID] = resid
258
259 // ZERO CROSSINGS, not first-difference extrema: crossings are robust to quantization where
260 // extrema are not. Amplitude for the decrement is the MAX EXCURSION WITHIN each half-cycle,
261 // which is the physically meaningful peak rather than one noisy sample.
262 var ncross: i64 = 0
263 var c1: i64 = 0 - 1
264 var c2: i64 = 0 - 1
265 var c3: i64 = 0 - 1
266 // PERIODIC-SPAN STATE. The span is grown while each successive crossing interval stays within
267 // one sample of the FIRST interval, which is one half period measured at the best signal-to-
268 // noise the record has. A crossing that is not a half-period marker -- the slow deadband graze
269 // AFTER the oscillation has died -- ends the span instead of poisoning it. That is the whole
270 // reason this estimator can REFUSE a corrupted record where the two-lobe form confidently
271 // publishes: MEASURED on an optical-flow series whose front end was injecting drift, the
272 // two-lobe decrement reported 78 permil against an injected 200 and this one abstained.
273 var iv0: i64 = 0 - 1
274 var span_ok: i64 = 0
275 var span_open: i64 = 1
276 var prev_c: i64 = 0 - 1
277 var lobe_pk: i64 = 0
278 var span_a: i64 = 0
279 var span_b: i64 = 0
280 var span_end: i64 = 0 - 1
281 var lobes: i64 = 0
282 var sgn: i64 = 0
283 var j: i64 = 0
284 while j < n {
285 var s2: i64 = 0
286 if sig[j] > quant { s2 = 1 }
287 if sig[j] < 0 - quant { s2 = 0 - 1 }
288 if s2 != 0 {
289 if sgn == 0 { sgn = s2 }
290 else {
291 if s2 != sgn {
292 ncross = ncross + 1
293 if c1 < 0 { c1 = j }
294 else { if c2 < 0 { c2 = j } else { if c3 < 0 { c3 = j } } }
295 sgn = s2
296 // ---- PERIODIC-SPAN TRACKING, in the SAME pass. Purely additive: every
297 // statement above is untouched and in the same order, so every field the
298 // incumbent published is computed exactly as it was.
299 if prev_c >= 0 {
300 lobes = lobes + 1
301 let iv: i64 = j - prev_c
302 if iv0 < 0 {
303 iv0 = iv
304 span_a = lobe_pk
305 span_b = lobe_pk
306 span_ok = 1
307 span_end = j
308 }
309 else {
310 if span_open == 1 {
311 if rd_abs(iv - iv0) <= RD_SPAN_TOL_SAMPLES {
312 span_ok = span_ok + 1
313 span_end = j
314 span_b = lobe_pk
315 }
316 else { span_open = 0 }
317 }
318 }
319 }
320 prev_c = j
321 lobe_pk = 0
322 }
323 }
324 }
325 if rd_abs(sig[j]) > lobe_pk { lobe_pk = rd_abs(sig[j]) }
326 j = j + 1
327 }
328 out[RD_R_PEAKS] = ncross
329 // ---- SPAN RESULTS, PUBLISHED BEFORE ANY EARLY RETURN BELOW, so a consumer always sees this
330 // estimator's OWN third state rather than a zero it cannot tell from a measurement.
331 out[RD_R_SPAN_LOBES] = lobes
332 out[RD_R_SPAN_RC] = RD_UNOBS_NO_PERIODIC_SPAN
333 if span_ok >= RD_SPAN_MIN_INTERVALS {
334 let halves: i64 = span_ok - 1
335 out[RD_R_SPAN_HALVES] = halves
336 out[RD_R_SPAN_A] = span_a
337 out[RD_R_SPAN_B] = span_b
338 let sp_us: i64 = RD_HALVES_PER_PERIOD * (span_end - c1) * dt_us / span_ok
339 out[RD_R_SPAN_PER_US] = sp_us
340 if sp_us > 0 { out[RD_R_SPAN_FN_MHZ] = RD_US_PER_S * RD_MHZ_PER_HZ / sp_us }
341 out[RD_R_SPAN_FN_TOL] = rd_span_fn_tol_permil(sp_us / dt_us, span_ok, span_b)
342 if span_a > span_b {
343 out[RD_R_SPAN_ZETA] = rd_zeta_permil_span(span_a, span_b, halves)
344 out[RD_R_SPAN_ZETA_TOL] = rd_span_zeta_tol_permil(span_a, span_b, halves)
345 out[RD_R_SPAN_RC] = RD_OK
346 }
347 else { out[RD_R_SPAN_RC] = RD_UNOBS_NO_DECAY }
348 }
349 if c1 <= 0 { return RD_UNOBS_NO_CROSSING }
350
351 var a1: i64 = 0
352 var a2: i64 = 0
353 var u: i64 = 0
354 while u < c1 { if rd_abs(sig[u]) > a1 { a1 = rd_abs(sig[u]) } u = u + 1 }
355 if c3 > 0 {
356 if c2 > 0 {
357 var v2: i64 = c2
358 while v2 < c3 { if rd_abs(sig[v2]) > a2 { a2 = rd_abs(sig[v2]) } v2 = v2 + 1 }
359 }
360 }
361 out[RD_R_A1] = a1
362 out[RD_R_A2] = a2
363 if c3 <= c1 { return RD_UNOBS_NO_PERIOD }
364
365 // c1 -> c3 spans exactly one full period (two sign changes)
366 let per_us: i64 = (c3 - c1) * dt_us
367 out[RD_R_PER_US] = per_us
368 if per_us > 0 { out[RD_R_FN_MHZ] = RD_US_PER_S * RD_MHZ_PER_HZ / per_us }
369 out[RD_R_ZETA] = rd_zeta_permil(a1, a2)
370 if a1 <= a2 { return RD_UNOBS_NO_DECAY }
371 return RD_OK
372}
373
374// ---- PHASE LAG between a DRIVE series and a RESPONSE series, by normalised cross-correlation
375// peak over a bounded lag search. Returned in PERMIL OF ONE DRIVE PERIOD, which is dimensionless
376// and therefore survives uncalibrated video where an absolute time lag would not mean anything
377// without knowing the subject's cadence. Returns RD_OK or a third-state reason.
378func rd_phase_lag_permil(drive: *i64, resp: *i64, n: i64, per_samples: i64, out_lag: *i64) -> i64 {
379 out_lag[0] = 0
380 if n <= 0 { return RD_UNOBS_NO_SAMPLES }
381 if per_samples <= 0 { return RD_UNOBS_NO_PERIOD }
382 var best: i64 = 0 - 1
383 var bestlag: i64 = 0
384 var found: i64 = 0
385 var lag: i64 = 0
386 while lag < per_samples {
387 var acc: i64 = 0
388 var cnt: i64 = 0
389 var i: i64 = 0
390 while i + lag < n {
391 acc = acc + drive[i] * resp[i + lag]
392 cnt = cnt + 1
393 i = i + 1
394 }
395 if cnt > 0 {
396 let norm: i64 = acc / cnt
397 if found == 0 { best = norm; bestlag = lag; found = 1 }
398 else { if norm > best { best = norm; bestlag = lag } }
399 }
400 lag = lag + 1
401 }
402 if found == 0 { return RD_UNOBS_NO_SAMPLES }
403 out_lag[0] = bestlag * RD_PERMIL / per_samples
404 return RD_OK
405}
406
407// ---- BASELINE REMOVAL FOR AN INTEGRATED SIGNAL -----------------------------------------------
408//
409// nx_dyna_bind_skin knows its equilibrium exactly: it settles the cage under gravity and subtracts
410// that level. A displacement series recovered from VIDEO has no settled state and is built by
411// INTEGRATING per-frame velocity, so a small per-frame bias accumulates into a ramp. Both removals
412// below RETURN what they removed, because a large drift or a large equilibrium shift is itself a
413// finding -- the region stopped tracking what it was pointed at -- and silently subtracting it
414// would convert a tracking failure into a clean-looking measurement.
415
416// FIT THE RAMP WHERE THE OSCILLATION HAS ALREADY SETTLED. Fitting over the WHOLE window is wrong
417// for a ringdown: a decaying cosine has a non-zero least-squares slope that is NOT drift, and
418// removing it tilts the waveform. MEASURED on a fixture whose true damped frequency is 1960 mHz:
419// a whole-window fit reported 7500 mHz because the tilt manufactured extra zero crossings near the
420// start. fit_from selects the settled region; the line is subtracted across the WHOLE series,
421// which is what makes it a drift removal rather than a windowed edit.
422func rd_detrend_from(sig: *i64, n: i64, fit_from: i64, out_a_q16: *i64, out_b_q16: *i64) -> i64 {
423 out_a_q16[0] = 0
424 out_b_q16[0] = 0
425 if n < 2 { return RD_UNOBS_NO_SAMPLES }
426 var f0: i64 = fit_from
427 if f0 < 0 { f0 = 0 }
428 if f0 > n - 2 { f0 = n - 2 }
429 var si: i64 = 0
430 var sx: i64 = 0
431 var six: i64 = 0
432 var si2: i64 = 0
433 var m: i64 = 0
434 var i: i64 = f0
435 while i < n {
436 si = si + i
437 sx = sx + sig[i]
438 six = six + i * sig[i]
439 si2 = si2 + i * i
440 m = m + 1
441 i = i + 1
442 }
443 if m < 2 { return RD_UNOBS_NO_SAMPLES }
444 let den: i64 = m * si2 - si * si
445 if den == 0 { return RD_UNOBS_NO_SAMPLES }
446 let b_q16: i64 = (m * six - si * sx) * RD_Q16 / den
447 let a_q16: i64 = (sx * RD_Q16 - b_q16 * si) / m
448 out_a_q16[0] = a_q16
449 out_b_q16[0] = b_q16
450 var j: i64 = 0
451 while j < n { sig[j] = sig[j] - (a_q16 + b_q16 * j) / RD_Q16; j = j + 1 }
452 return RD_OK
453}
454
455// Whole-window fit: the ordinary contract, expressed over the scoped one.
456func rd_detrend(sig: *i64, n: i64, out_a_q16: *i64, out_b_q16: *i64) -> i64 {
457 return rd_detrend_from(sig, n, 0, out_a_q16, out_b_q16)
458}
459
460// A free decay settles at a NEW equilibrium, not at zero, and a LINEAR fit cannot remove a STEP.
461// The estimator is the mean of the final eighth -- the same window rd_analyze already treats as the
462// residual, so there is ONE definition of "the tail" in this file rather than two.
463func rd_center_on_tail(sig: *i64, n: i64, out_eq: *i64) -> i64 {
464 out_eq[0] = 0
465 if n < 8 { return RD_UNOBS_NO_SAMPLES }
466 var acc: i64 = 0
467 var cnt: i64 = 0
468 var t: i64 = n - n/8
469 while t < n { acc = acc + sig[t]; cnt = cnt + 1; t = t + 1 }
470 if cnt <= 0 { return RD_UNOBS_NO_SAMPLES }
471 let eq: i64 = acc / cnt
472 out_eq[0] = eq
473 var j: i64 = 0
474 while j < n { sig[j] = sig[j] - eq; j = j + 1 }
475 return RD_OK
476}