|
| 1 | +# Minimal resonance + phi-lock helpers for TFT |
| 2 | +# Deterministic, numpy-first; scipy used if present. |
| 3 | + |
| 4 | +from __future__ import annotations |
| 5 | +import numpy as np |
| 6 | + |
| 7 | +def rng(seed: int | None = 1337) -> np.random.Generator: |
| 8 | + return np.random.default_rng(seed) |
| 9 | + |
| 10 | +def rotation_matrix(dim: int = 3, seed: int | None = 1337) -> np.ndarray: |
| 11 | + """Random orthogonal R (det=+1) via QR, seeded.""" |
| 12 | + g = rng(seed) |
| 13 | + A = g.normal(size=(dim, dim)) |
| 14 | + Q, _ = np.linalg.qr(A) |
| 15 | + if np.linalg.det(Q) < 0: |
| 16 | + Q[:, 0] *= -1 |
| 17 | + return Q |
| 18 | + |
| 19 | +def spd_matrix(dim: int = 3, seed: int | None = 1337) -> np.ndarray: |
| 20 | + """Symmetric positive-definite matrix A A^T, seeded.""" |
| 21 | + g = rng(seed) |
| 22 | + A = g.normal(size=(dim, dim)) |
| 23 | + return A @ A.T |
| 24 | + |
| 25 | +def rotate_rank2(T: np.ndarray, R: np.ndarray) -> np.ndarray: |
| 26 | + """Orthogonal congruence: T' = R T R^T (preserves eigvals, Fro norm).""" |
| 27 | + return R @ T @ R.T |
| 28 | + |
| 29 | +def invariants(T: np.ndarray) -> dict: |
| 30 | + """Numerical invariants we care about for the demo.""" |
| 31 | + # Symmetrize to be safe numerically |
| 32 | + S = 0.5 * (T + T.T) |
| 33 | + fro = np.linalg.norm(S, "fro") |
| 34 | + # eigvalsh: Hermitian (real symmetric) eigenvalues |
| 35 | + ev = np.linalg.eigvalsh(S) |
| 36 | + return {"fro": float(fro), "eigvals": np.sort(ev)} |
| 37 | + |
| 38 | +def _analytic_signal_fft(x: np.ndarray, axis: int = -1) -> np.ndarray: |
| 39 | + """ |
| 40 | + Hilbert analytic signal via FFT (no scipy required). |
| 41 | + Produces y = x + i*x_hilbert with 90° quadrature. |
| 42 | + """ |
| 43 | + X = np.fft.fft(x, axis=axis) |
| 44 | + N = x.shape[axis] |
| 45 | + H = np.zeros(N) |
| 46 | + if N % 2 == 0: |
| 47 | + H[0] = 1 |
| 48 | + H[N//2] = 1 |
| 49 | + H[1:N//2] = 2 |
| 50 | + else: |
| 51 | + H[0] = 1 |
| 52 | + H[1:(N+1)//2] = 2 |
| 53 | + shape = [1] * x.ndim |
| 54 | + shape[axis] = N |
| 55 | + H = H.reshape(shape) |
| 56 | + return np.fft.ifft(X * H, axis=axis) |
| 57 | + |
| 58 | +def analytic_signal(x: np.ndarray, axis: int = -1) -> np.ndarray: |
| 59 | + """Try scipy.signal.hilbert; fall back to FFT implementation.""" |
| 60 | + try: |
| 61 | + from scipy.signal import hilbert # type: ignore |
| 62 | + return hilbert(x, axis=axis) |
| 63 | + except Exception: |
| 64 | + return _analytic_signal_fft(x, axis=axis) |
| 65 | + |
| 66 | +def phi_lock_pair(x: np.ndarray, axis: int = -1) -> tuple[np.ndarray, np.ndarray]: |
| 67 | + """ |
| 68 | + Return (real, imag) quadrature pair at 90° (φ = π/2). |
| 69 | + Energies normalized to match. |
| 70 | + """ |
| 71 | + z = analytic_signal(x, axis=axis) |
| 72 | + a = np.real(z) |
| 73 | + b = np.imag(z) |
| 74 | + # Normalize energies to match (avoid trivial scale mismatch) |
| 75 | + ea = np.sqrt(np.mean(a**2)) |
| 76 | + eb = np.sqrt(np.mean(b**2)) + 1e-12 |
| 77 | + b *= (ea / eb) |
| 78 | + return a, b |
| 79 | + |
| 80 | +def map_to_audio(eigvals: np.ndarray, fmin=220.0, fmax=880.0) -> np.ndarray: |
| 81 | + """ |
| 82 | + Map sorted eigenvalues to audible frequencies, linearly. |
| 83 | + Handles degenerate cases gracefully. |
| 84 | + """ |
| 85 | + v = np.real(eigvals) |
| 86 | + vmin, vmax = float(np.min(v)), float(np.max(v)) |
| 87 | + if np.isclose(vmax, vmin): |
| 88 | + return np.full_like(v, (fmin + fmax) * 0.5, dtype=float) |
| 89 | + return fmin + (v - vmin) * (fmax - fmin) / (vmax - vmin) |
0 commit comments