#!/usr/bin/env python3 """ OVERTONE ANALYSIS — Dimensional Boundary Harmonic Detection ================================================================================ Tests the overtone hypothesis: do dimensional boundaries produce harmonic series in the temporal pulse sequence? METHOD: 1. Run v6c simulation and record EVERY pulse event (time, energy, dimension) 2. Compute the inter-pulse interval (IPI) series for each cascade transition 3. FFT the IPI series to find frequency components 4. Check: does the temporal pattern contain harmonics at 2x, 3x, 5x, 7x, 11x the fundamental injection frequency? 5. Measure harmonic amplitude vs harmonic number — does it follow the expected 1/n attenuation of an overtone series? PREDICTION (overtone hypothesis): - The pulse sequence should show harmonic content beyond the fundamental - Harmonics at {2,3} should be strongest - {5} should be present but weaker - {7} should be present but weaker still - {11} should be faint - Amplitude should decrease with harmonic number NULL HYPOTHESIS: - The pulse sequence is periodic at one frequency (the overflow frequency) - No significant harmonic content beyond the fundamental - The dimensional boundary acts as a simple frequency converter, not an overtone generator If the null holds, the overtone hypothesis loses its computational support. """ import json import math import numpy as np import time as time_mod from datetime import datetime, timezone from pathlib import Path import sys sys.path.insert(0, str(Path(__file__).parent)) from importlib import import_module spec = import_module("SIM-003_v6c_wave_reinjection") CascadeEngine = spec.CascadeEngine RESULTS_DIR = Path(__file__).parent / "results_overtone_analysis" RESULTS_DIR.mkdir(exist_ok=True) def run_with_pulse_recording(element="Fe", n_cycles=3000, N3=128): """Run v6c and record every pulse event with precise timing.""" engine = CascadeEngine(element, delta_r_base=0.0001, N1=2048, N2=256, N3=N3, max_dim=3) spc = int(engine.T_pulse / engine.dt_master) total_steps = n_cycles * spc # Pulse event logs: (step, time, energy) pulses_1to2 = [] # 1D overflow → 2D injection pulses_2to3 = [] # 2D overflow → 3D injection eq_measured = {1: False, 2: False, 3: False} print(f"Running {element}, {n_cycles} cycles, {N3}³ grid...") print(f"Recording all pulse events for overtone analysis.") t_start = time_mod.monotonic() for step in range(total_steps): t = step * engine.dt_master # 1D overflow_1d = engine.grids[1].step(t) if 1 in engine.grids else 0 if not eq_measured[1] and engine.grids[1].is_equilibrium: engine._measure_overflow_quantum(1) eq_measured[1] = True if overflow_1d > 0: engine.overflow_buffer[1] += overflow_1d # 2D overflow_2d = 0.0 if 2 in engine.grids: q1 = engine.overflow_quantum[1] if q1 is not None: if not engine.grids[2].active and engine.overflow_buffer[1] >= q1: engine.grids[2].active = True if engine.grids[2].active: if not engine.grids[2].injecting and engine.overflow_buffer[1] >= q1: engine.grids[2].start_cascade_pulse(q1, t) engine.overflow_buffer[1] -= q1 engine.pulse_count[1] += 1 # RECORD the pulse event pulses_1to2.append((step, t, float(q1))) for _ in range(max(1, int(engine.grids[2].dt / engine.dt_master))): overflow_2d = engine.grids[2].step(t) if not eq_measured[2] and engine.grids[2].is_equilibrium: engine._measure_overflow_quantum(2) eq_measured[2] = True if overflow_2d > 0: engine.overflow_buffer[2] += overflow_2d # 3D if 3 in engine.grids: q2 = engine.overflow_quantum[2] if q2 is not None: if not engine.grids[3].active and engine.overflow_buffer[2] >= q2: engine.grids[3].active = True if engine.grids[3].active: if not engine.grids[3].injecting and engine.overflow_buffer[2] >= q2: engine.grids[3].start_cascade_pulse(q2, t) engine.overflow_buffer[2] -= q2 engine.pulse_count[2] += 1 # RECORD the pulse event pulses_2to3.append((step, t, float(q2))) for _ in range(max(1, int(engine.grids[3].dt / engine.dt_master))): engine.grids[3].step(t) if not eq_measured[3] and engine.grids[3].is_equilibrium: engine._measure_overflow_quantum(3) eq_measured[3] = True # Progress if step > 0 and step % (total_steps // 10) == 0: cycle = step / spc print(f" c={cycle:.0f} pulses: 1->2={len(pulses_1to2)} 2->3={len(pulses_2to3)}") elapsed = time_mod.monotonic() - t_start print(f"Done in {elapsed:.1f}s. Pulses: 1->2={len(pulses_1to2)}, 2->3={len(pulses_2to3)}") return { "element": element, "n_cycles": n_cycles, "N3": N3, "elapsed_s": elapsed, "dt_master": engine.dt_master, "pulses_1to2": pulses_1to2, "pulses_2to3": pulses_2to3, "final_state": { "1D_global_r": engine.grids[1].global_r, "2D_global_r": engine.grids[2].global_r if 2 in engine.grids else 0, "3D_global_r": engine.grids[3].global_r if 3 in engine.grids else 0, } } def analyze_overtones(pulse_events, dt_master, label=""): """Analyze a pulse sequence for harmonic content. pulse_events: list of (step, time, energy) tuples Returns: overtone analysis results """ if len(pulse_events) < 10: return {"label": label, "error": "Too few pulses", "n_pulses": len(pulse_events)} times = np.array([p[1] for p in pulse_events]) energies = np.array([p[2] for p in pulse_events]) # Inter-pulse intervals ipis = np.diff(times) mean_ipi = np.mean(ipis) std_ipi = np.std(ipis) cv_ipi = std_ipi / mean_ipi if mean_ipi > 0 else 0 # Fundamental frequency from mean IPI f_fundamental = 1.0 / mean_ipi if mean_ipi > 0 else 0 print(f"\n=== OVERTONE ANALYSIS: {label} ===") print(f" Pulses: {len(pulse_events)}") print(f" Mean IPI: {mean_ipi:.6f}") print(f" Std IPI: {std_ipi:.6f}") print(f" CV: {cv_ipi:.4f}") print(f" Fundamental freq: {f_fundamental:.4f}") # FFT of the inter-pulse interval series # This reveals periodic modulations of the pulse timing if len(ipis) < 20: return {"label": label, "error": "Too few intervals for FFT", "n_pulses": len(pulse_events), "mean_ipi": mean_ipi} # Subtract mean to focus on variations ipi_centered = ipis - mean_ipi fft = np.abs(np.fft.rfft(ipi_centered)) freqs = np.fft.rfftfreq(len(ipi_centered)) # Also FFT the energy sequence — amplitude modulation energy_centered = energies[1:] - np.mean(energies[1:]) # skip first, align with IPI if len(energy_centered) == len(ipi_centered): fft_energy = np.abs(np.fft.rfft(energy_centered)) else: fft_energy = np.zeros_like(fft) # Normalize FFTs fft_norm = fft / np.max(fft) if np.max(fft) > 0 else fft fft_e_norm = fft_energy / np.max(fft_energy) if np.max(fft_energy) > 0 else fft_energy # Look for peaks at harmonic ratios of the fundamental # In the IPI FFT, frequency is in units of "cycles per pulse" # Harmonic n of the fundamental appears at freq = n / (total_pulses) # But more useful: look at the actual frequency peaks # Find significant peaks threshold = 0.05 peaks_ipi = [] for i in range(1, len(fft_norm) - 1): if (fft_norm[i] > fft_norm[i-1] and fft_norm[i] > fft_norm[i+1] and fft_norm[i] > threshold): peaks_ipi.append({ "index": int(i), "freq": float(freqs[i]), "amplitude": float(fft[i]), "normalized": float(fft_norm[i]), }) peaks_ipi.sort(key=lambda p: -p["normalized"]) peaks_energy = [] for i in range(1, len(fft_e_norm) - 1): if (fft_e_norm[i] > fft_e_norm[i-1] and fft_e_norm[i] > fft_e_norm[i+1] and fft_e_norm[i] > threshold): peaks_energy.append({ "index": int(i), "freq": float(freqs[i]), "amplitude": float(fft_energy[i]), "normalized": float(fft_e_norm[i]), }) peaks_energy.sort(key=lambda p: -p["normalized"]) print(f"\n IPI FFT peaks (>{threshold:.0%} of max):") if peaks_ipi: for p in peaks_ipi[:10]: # Express frequency as ratio to the lowest peak print(f" freq={p['freq']:.4f} amp={p['normalized']:.3f}") else: print(f" NONE — pulse timing is uniform (supports null hypothesis)") print(f"\n Energy FFT peaks:") if peaks_energy: for p in peaks_energy[:10]: print(f" freq={p['freq']:.4f} amp={p['normalized']:.3f}") else: print(f" NONE — pulse energy is constant") # Check for harmonic ratios between peaks harmonic_ratios = [] if len(peaks_ipi) >= 2: f0 = peaks_ipi[0]["freq"] if f0 > 0: print(f"\n Harmonic ratios (relative to strongest peak f={f0:.4f}):") for p in peaks_ipi[1:10]: ratio = p["freq"] / f0 nearest_int = round(ratio) deviation = abs(ratio - nearest_int) / nearest_int if nearest_int > 0 else 1 is_harmonic = deviation < 0.1 # within 10% of integer ratio marker = " <<< HARMONIC" if is_harmonic else "" harmonic_ratios.append({ "freq": p["freq"], "ratio": float(ratio), "nearest_integer": int(nearest_int), "deviation": float(deviation), "is_harmonic": is_harmonic, "amplitude": p["normalized"], }) print(f" f={p['freq']:.4f} ratio={ratio:.2f} " f"(nearest int: {nearest_int}, dev: {deviation:.2%}) " f"amp={p['normalized']:.3f}{marker}") # Check specific harmonic numbers: are {2,3,5,7,11} present? print(f"\n Checking for specific harmonics of fundamental:") target_harmonics = [2, 3, 5, 7, 11] harmonic_amplitudes = {} if len(peaks_ipi) > 0: f0 = peaks_ipi[0]["freq"] for n in target_harmonics: target_freq = f0 * n # Find nearest FFT bin if target_freq < freqs[-1]: idx = np.argmin(np.abs(freqs - target_freq)) amp = float(fft_norm[idx]) harmonic_amplitudes[n] = amp present = "PRESENT" if amp > 0.05 else "absent" print(f" {n}x fundamental (f={target_freq:.4f}): " f"amp={amp:.4f} [{present}]") else: harmonic_amplitudes[n] = 0 print(f" {n}x fundamental: beyond Nyquist") # Overall assessment n_harmonics_found = sum(1 for a in harmonic_amplitudes.values() if a > 0.05) print(f"\n ASSESSMENT:") if n_harmonics_found == 0 and len(peaks_ipi) == 0: print(f" NULL HYPOTHESIS SUPPORTED — no harmonic content") print(f" Pulse timing is uniform. Boundary does not produce overtones.") assessment = "null_supported" elif n_harmonics_found == 0 and len(peaks_ipi) > 0: print(f" PEAKS PRESENT but not at harmonic ratios") print(f" Some temporal structure, but not an overtone series.") assessment = "structure_non_harmonic" elif n_harmonics_found <= 2: print(f" WEAK HARMONIC CONTENT — {n_harmonics_found} harmonics detected") print(f" Suggestive but not definitive.") assessment = "weak_harmonic" else: print(f" HARMONIC SERIES DETECTED — {n_harmonics_found} of {len(target_harmonics)} harmonics present") print(f" Dimensional boundary produces overtones.") assessment = "harmonics_detected" return { "label": label, "n_pulses": len(pulse_events), "mean_ipi": float(mean_ipi), "std_ipi": float(std_ipi), "cv_ipi": float(cv_ipi), "f_fundamental": float(f_fundamental), "peaks_ipi": peaks_ipi[:20], "peaks_energy": peaks_energy[:20], "harmonic_ratios": harmonic_ratios, "harmonic_amplitudes": harmonic_amplitudes, "n_harmonics_found": n_harmonics_found, "assessment": assessment, "ipi_fft": [float(x) for x in fft_norm[:100]], "freqs": [float(x) for x in freqs[:100]], } def main(): import argparse parser = argparse.ArgumentParser(description="Overtone Analysis") parser.add_argument("--element", default="Fe") parser.add_argument("--cycles", type=int, default=3000) parser.add_argument("--grid-3d", type=int, default=128) args = parser.parse_args() print("=" * 70) print("OVERTONE ANALYSIS — DIMENSIONAL BOUNDARY HARMONICS") print(f" {datetime.now(timezone.utc).isoformat()}") print(f" Testing: do dimensional boundaries produce harmonic series?") print(f" Element: {args.element}, Cycles: {args.cycles}, Grid: {args.grid_3d}³") print("=" * 70) # Run simulation with pulse recording data = run_with_pulse_recording(args.element, args.cycles, args.grid_3d) # Analyze overtones in each cascade transition results = { "timestamp": datetime.now(timezone.utc).isoformat(), "element": args.element, "n_cycles": args.cycles, "grid_3d": args.grid_3d, "final_state": data["final_state"], } results["1D_to_2D"] = analyze_overtones( data["pulses_1to2"], data["dt_master"], label="1D→2D cascade") results["2D_to_3D"] = analyze_overtones( data["pulses_2to3"], data["dt_master"], label="2D→3D cascade") # Summary print(f"\n{'='*70}") print("SUMMARY") print(f"{'='*70}") for key in ["1D_to_2D", "2D_to_3D"]: r = results[key] print(f" {r['label']}: {r.get('assessment', 'N/A')} " f"({r.get('n_harmonics_found', 0)} harmonics, " f"{r.get('n_pulses', 0)} pulses)") # Save ts = datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S") out = RESULTS_DIR / f"overtone_{args.element}_{ts}.json" with open(out, "w") as f: json.dump(results, f, indent=2, default=float) print(f"\nResults: {out}") if __name__ == "__main__": main()