Create breath_keeper.py

breath_keeper.py

@@ -0,0 +1,219 @@
+"""
+Breath Keeper A/B analysis and persistence metrics.
+
+This module provides utilities for analyzing oscillator signals, including
+calculation of unbiased autocorrelation, coherence half-life, beat period,
+phase-slip, and energy drift. It also implements a breath keeper (phase-locked
+loop with node snapping, amplitude control, and symplectic oscillator) to
+maintain phase coherence and conserve energy. A command-line interface is
+provided for running baseline vs keeper-enabled analysis on a CSV file of
+time series data.
+"""
+
+from __future__ import annotations
+import numpy as np
+from dataclasses import dataclass
+from typing import Tuple, Optional, Dict
+
+try:
+    from scipy.signal import hilbert, find_peaks
+except Exception as e:
+    raise SystemExit("Please install scipy via 'pip install scipy' to use breath_keeper.") from e
+
+# Utility functions
+def unbiased_autocorr(x: np.ndarray) -> np.ndarray:
+    x = x.astype(np.float64)
+    x = x - np.mean(x)
+    N = len(x)
+    ac = np.correlate(x, x, mode='full')
+    ac = ac[N-1:]
+    norm = np.arange(N, 0, -1, dtype=np.float64)
+    ac_unbiased = ac / norm
+    return ac_unbiased / ac_unbiased[0]
+
+def fit_coherence_half_life(ac: np.ndarray, Fs: float, min_lag_s: float = 2.0, max_lag_frac: float = 0.5) -> Tuple[float, Tuple[float, float]]:
+    N = len(ac)
+    t = np.arange(N)/Fs
+    lo = int(min_lag_s*Fs)
+    hi = int(min(N-1, max_lag_frac*N))
+    if hi <= lo+5:
+        return float('nan'), (float('nan'), float('nan'))
+    y = np.clip(ac[lo:hi], 1e-12, 1.0)
+    tt = t[lo:hi]
+    A = np.vstack([tt, np.ones_like(tt)]).T
+    coeff, _, _, _ = np.linalg.lstsq(A, np.log(y), rcond=None)
+    slope, intercept = coeff
+    tau = -1.0 / slope if slope < 0 else float('nan')
+    residuals = np.log(y) - (A @ coeff)
+    sigma = np.std(residuals)
+    denom = np.sum((tt - np.mean(tt))**2)
+    if denom <= 0:
+        return tau, (float('nan'), float('nan'))
+    var_slope = sigma**2 / denom
+    se_slope = np.sqrt(var_slope)
+    slope_lo = slope - 2*se_slope
+    slope_hi = slope + 2*se_slope
+    tau_lo = -1.0/slope_hi if slope_hi < 0 else float('nan')
+    tau_hi = -1.0/slope_lo if slope_lo < 0 else float('nan')
+    return float(tau), (float(tau_lo), float(tau_hi))
+
+def analytic_envelope(x: np.ndarray) -> np.ndarray:
+    return np.abs(hilbert(x))
+
+def beat_period_from_envelope(env: np.ndarray, Fs: float, min_period_s: float = 0.5, max_period_s: float = 10.0) -> float:
+    ac = unbiased_autocorr(env)
+    lags = np.arange(len(ac))/Fs
+    lo = int(min_period_s*Fs)
+    hi = int(min(len(ac)-1, max_period_s*Fs))
+    if hi <= lo+3:
+        return float('nan')
+    peaks, _ = find_peaks(ac[lo:hi], height=0.2)
+    if len(peaks) == 0:
+        return float('nan')
+    first = peaks[0] + lo
+    return lags[first]
+
+def node_trough_indices(env: np.ndarray, Tb: float, Fs: float) -> np.ndarray:
+    if not np.isfinite(Tb) or Tb <= 0:
+        idx, _ = find_peaks(-env, distance=int(0.25*Fs))
+        return idx
+    idx, _ = find_peaks(-env, distance=int(0.8*Tb*Fs))
+    return idx
+
+def phase_from_analytic(x: np.ndarray) -> np.ndarray:
+    return np.unwrap(np.angle(hilbert(x)))
+
+def energy_series(x: np.ndarray, Fs: float, win_periods: float = 1.0, carrier_Hz: Optional[float] = None) -> np.ndarray:
+    if carrier_Hz and carrier_Hz > 0:
+        win = int(max(1, round(Fs/carrier_Hz*win_periods)))
+    else:
+        win = int(max(1, round(Fs/10)))
+    kernel = np.ones(win)/win
+    rms = np.sqrt(np.convolve(x**2, kernel, mode='same'))
+    return rms**2
+
+def wrap_pi(a):
+    return (a + np.pi) % (2*np.pi) - np.pi
+
+@dataclass
+class KeeperConfig:
+    Fs: float
+    Kp: float = 0.05
+    Ki: float = 0.001
+    agc_tau: float = 0.5
+    snap_thresh: float = 0.05
+    omega_init: Optional[float] = None
+
+class BreathKeeper:
+    def __init__(self, cfg: KeeperConfig):
+        self.cfg = cfg
+        self._phi_int = 0.0
+        self._amp = 1.0
+        self.q = 0.0
+        self.p = 0.0
+        self.omega = cfg.omega_init if cfg.omega_init else 2*np.pi*1.0
+
+    def phase_est(self) -> float:
+        return np.arctan2(self.p, self.q + 1e-12)
+
+    def set_phase(self, phi: float):
+        A = max(1e-9, self._amp)
+        self.q = A * np.cos(phi)
+        self.p = A * np.sin(phi)
+
+    def step(self, x_t: float, env_t: float, phi_meas: float) -> float:
+        phi_err = wrap_pi(phi_meas - self.phase_est())
+        self._phi_int += self.cfg.Ki * phi_err
+        dphi = self.cfg.Kp * phi_err + self._phi_int
+        self.omega = max(1e-6, self.omega + dphi)
+        if env_t < self.cfg.snap_thresh * (self._amp + 1e-9):
+            self.set_phase(np.round(self.phase_est()/np.pi)*np.pi)
+        alpha = np.exp(-1.0/(self.cfg.agc_tau*self.cfg.Fs))
+        self._amp = alpha*self._amp + (1-alpha)*abs(x_t)
+        dt = 1.0/self.cfg.Fs
+        self.p -= (self.omega**2)*self.q*(dt*0.5)
+        self.q += self.p*dt
+        self.p -= (self.omega**2)*self.q*(dt*0.5)
+        return self.phase_est()
+
+@dataclass
+class Metrics:
+    Tb: float
+    tau_c: float
+    tau_ci: Tuple[float,float]
+    phase_slip_rad_per_beat: float
+    energy_drift_per_beat: float
+
+def compute_metrics(x: np.ndarray, Fs: float) -> Metrics:
+    env = analytic_envelope(x)
+    Tb = beat_period_from_envelope(env, Fs)
+    ac = unbiased_autocorr(env)
+    tau_c, tau_ci = fit_coherence_half_life(ac, Fs)
+    nodes = node_trough_indices(env, Tb, Fs)
+    phi = phase_from_analytic(x)
+    if len(nodes) >= 2:
+        dphi = wrap_pi(np.diff(phi[nodes]))
+        phase_slip = float(np.mean(np.abs(dphi)))
+    else:
+        phase_slip = float('nan')
+    E = energy_series(x, Fs)
+    if len(nodes) >= 2:
+        Eb = E[nodes]
+        rel_changes = np.diff(Eb)/(Eb[:-1]+1e-12)
+        energy_drift = float(np.mean(rel_changes))
+    else:
+        energy_drift = float('nan')
+    return Metrics(Tb=Tb, tau_c=float(tau_c), tau_ci=(float(tau_ci[0]), float(tau_ci[1])), phase_slip_rad_per_beat=phase_slip, energy_drift_per_beat=energy_drift)
+
+def keeper_follow(x: np.ndarray, Fs: float) -> np.ndarray:
+    env = analytic_envelope(x)
+    phi_meas = np.unwrap(np.angle(hilbert(x)))
+    cfg = KeeperConfig(Fs=Fs)
+    k = BreathKeeper(cfg)
+    y = np.zeros_like(x)
+    for n in range(len(x)):
+        phi = k.step(x[n], env[n], phi_meas[n])
+        y[n] = np.cos(phi)
+    return y
+
+def analyze_ab(x: np.ndarray, Fs: float) -> Dict[str, Metrics]:
+    base_metrics = compute_metrics(x, Fs)
+    y = keeper_follow(x, Fs)
+    keep_metrics = compute_metrics(y, Fs)
+    return {"baseline": base_metrics, "keeper": keep_metrics}
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(description="Breath keeper analysis: compute baseline and keeper metrics on CSV file.")
+    parser.add_argument("--csv", type=str, required=True, help="Path to CSV with columns t,x or x.")
+    parser.add_argument("--fs", type=float, required=False, help="Sampling rate in Hz if no t column.")
+    args = parser.parse_args()
+    import numpy as np
+    data = np.genfromtxt(args.csv, delimiter=",", names=True, dtype=None, encoding=None)
+    if "x" in data.dtype.names:
+        x = data["x"].astype(np.float64)
+        if "t" in data.dtype.names:
+            t = data["t"].astype(np.float64)
+            Fs_val = 1.0/np.mean(np.diff(t))
+        else:
+            if args.fs is None:
+                raise ValueError("Sampling rate must be provided if no t column.")
+            Fs_val = float(args.fs)
+    else:
+        raw = np.genfromtxt(args.csv, delimiter=",")
+        if raw.ndim == 1:
+            if args.fs is None:
+                raise ValueError("Sampling rate must be provided for single column CSV.")
+            x = raw.astype(np.float64)
+            Fs_val = float(args.fs)
+        else:
+            t = raw[:,0].astype(np.float64)
+            x = raw[:,1].astype(np.float64)
+            Fs_val = 1.0/np.mean(np.diff(t))
+    results = analyze_ab(x, Fs_val)
+    for label, m in results.items():
+        print(f"[{label}]")
+        print(f"  Beat period Tb (s):         {m.Tb:.6f}")
+        print(f"  Coherence half-life \u03c4c (s): {m.tau_c:.6f}  (CI ~ {m.tau_ci[0]:.3f}, {m.tau_ci[1]:.3f})")
+        print(f"  Phase slip |\u03c6̇| (rad/beat):  {m.phase_slip_rad_per_beat:.6e}")
+        print(f"  Energy drift \u0110 (/beat):      {m.energy_drift_per_beat:.6e}")