# %% [markdown]
# # FRCFD Deep Dive – Single Cell (Light Version)
#
# This cell:
# - Installs `gwpy` if needed
# - Fetches 0.5 seconds of H1 and L1 data before GW150914
# - Extracts the 200.2 Hz amplitude envelope via heterodyne demodulation
# - Computes cross‑correlation between H1 and L1 envelopes
# - Plots envelopes vs. expected orbital frequency
# - Performs phase‑amplitude coupling analysis
# - Computes FFT of the envelope to look for orbital modulation
# - Prints a comprehensive summary of findings
# ---------------------------
# 0. Setup and install
# ---------------------------
!pip install gwpy -q
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt, hilbert, correlate
from gwpy.timeseries import TimeSeries
import gc
import warnings
warnings.filterwarnings("ignore")
print("Imports done.\n")
# ---------------------------
# 1. Parameters (light mode)
# ---------------------------
M_CHIRP = 28.1 # solar masses (GW150914 final estimate)
G = 6.67430e-11
c = 299792458
M_sun = 1.989e30
M_chirp_kg = M_CHIRP * M_sun
F_ANCHOR = 200.2 # Hz (FRCFD invariant)
BANDWIDTH = 10 # Hz (195–205 Hz)
DURATION = 0.5 # seconds – shorter to save memory
T0 = 1126259462.423 # merger peak GPS
START_TIME = T0 - DURATION
END_TIME = T0 - 0.05 # stop 50 ms before peak
DELAY_H1_L1 = 0.0069 # seconds (light travel time)
print(f"Parameters set: {DURATION}s data, {F_ANCHOR} Hz anchor, chirp mass {M_CHIRP} Msun\n")
# ---------------------------
# 2. Orbital frequency function
# ---------------------------
def orbital_frequency(t_rel):
tau = -t_rel
tau = np.maximum(tau, 1e-6)
GMc3 = (G * M_chirp_kg) / c**3
factor = (5 / 256) * (1 / tau)
f_orb = (1 / np.pi) * (factor)**(3/8) * (GMc3)**(-5/8)
return f_orb
# ---------------------------
# 3. Helper: extract envelope for one detector
# ---------------------------
def extract_envelope_light(detector, start, end, f_center, bw):
data = TimeSeries.fetch_open_data(detector, start, end, verbose=False)
fs = data.sample_rate.value
strain = data.value
t = np.arange(len(strain)) / fs + start - T0 # time relative to merger
# Bandpass filter
nyq = fs / 2
low = max(f_center - bw/2, 0.1)
high = f_center + bw/2
b, a = butter(4, [low/nyq, high/nyq], btype='band')
filtered = filtfilt(b, a, strain)
# Heterodyne demodulation
t_vec = np.arange(len(filtered)) / fs
analytic = filtered * np.exp(-1j * 2 * np.pi * f_center * t_vec)
b_lp, a_lp = butter(4, 200/nyq, btype='low')
env_complex = filtfilt(b_lp, a_lp, analytic)
envelope = np.abs(env_complex)
# Clean up
del data, strain, filtered, analytic, env_complex
gc.collect()
return t, envelope, fs
# ---------------------------
# 4. Extract H1 and L1 envelopes
# ---------------------------
print("Extracting H1 envelope...")
t_h1, env_h1, fs = extract_envelope_light('H1', START_TIME, END_TIME, F_ANCHOR, BANDWIDTH)
print("Extracting L1 envelope...")
t_l1, env_l1, fs_l1 = extract_envelope_light('L1', START_TIME, END_TIME, F_ANCHOR, BANDWIDTH)
# Ensure same time axis
if len(t_h1) != len(t_l1):
min_len = min(len(t_h1), len(t_l1))
t_h1 = t_h1[:min_len]
env_h1 = env_h1[:min_len]
t_l1 = t_l1[:min_len]
env_l1 = env_l1[:min_len]
t = t_h1 # use common time
print(f"Data length: {len(t)} samples, fs = {fs:.1f} Hz\n")
# ---------------------------
# 5. Cross‑correlation (H1 vs L1)
# ---------------------------
env_h1_ac = env_h1 - np.mean(env_h1)
env_l1_ac = env_l1 - np.mean(env_l1)
corr = correlate(env_h1_ac, env_l1_ac, mode='same')
lags = np.linspace(-len(env_h1_ac)//2, len(env_h1_ac)//2, len(corr)) / fs
plt.figure(figsize=(10,4))
plt.plot(lags*1000, corr)
plt.axvline(DELAY_H1_L1*1000, color='r', linestyle='--', label=f'Expected {DELAY_H1_L1*1000:.1f} ms')
plt.axvline(-DELAY_H1_L1*1000, color='r', linestyle='--')
plt.xlabel('Lag (ms)')
plt.ylabel('Cross-correlation')
plt.title('H1 vs L1 – 200.2 Hz Envelope Cross-correlation')
plt.legend()
plt.grid(True)
plt.show()
peak_idx = np.argmax(np.abs(corr))
peak_lag_ms = lags[peak_idx] * 1000
print(f"Peak cross-correlation at lag = {peak_lag_ms:.2f} ms")
lag_match = np.abs(peak_lag_ms - DELAY_H1_L1*1000) < 5
print(f"Matches expected {DELAY_H1_L1*1000:.1f} ms? {'YES' if lag_match else 'NO'}\n")
# ---------------------------
# 6. Plot envelopes + orbital frequency
# ---------------------------
f_orb = orbital_frequency(t)
plt.figure(figsize=(12,5))
plt.plot(t, env_h1, label='H1 envelope', alpha=0.7)
plt.plot(t, env_l1, label='L1 envelope', alpha=0.7)
plt.xlabel('Time relative to merger (s)')
plt.ylabel('200.2 Hz Amplitude')
plt.legend()
plt.grid(True)
ax2 = plt.twinx()
ax2.plot(t, f_orb, 'r--', label='Orbital freq (chirp)')
ax2.set_ylabel('Frequency (Hz)', color='red')
plt.title('200.2 Hz Envelopes vs. Expected Orbital Frequency')
plt.show()
# ---------------------------
# 7. Phase‑Amplitude Coupling (PAC)
# ---------------------------
# Fetch broadband H1 data for phase
print("Fetching broadband H1 data for phase...")
data_gw = TimeSeries.fetch_open_data('H1', START_TIME, END_TIME, verbose=False)
strain_gw = data_gw.value
b_broad, a_broad = butter(4, [30/fs, 400/fs], btype='band')
gw_broad = filtfilt(b_broad, a_broad, strain_gw)
phase_gw = np.angle(hilbert(gw_broad))
del data_gw, strain_gw, gw_broad
gc.collect()
# Align phase to envelope time axis
if len(phase_gw) != len(t):
phase_gw = np.interp(t, np.arange(len(phase_gw))/fs + START_TIME - T0, phase_gw)
# Bin envelope amplitude by phase
n_bins = 18
bins = np.linspace(-np.pi, np.pi, n_bins+1)
bin_centers = (bins[:-1] + bins[1:])/2
mean_amp = []
for i in range(n_bins):
mask = (phase_gw >= bins[i]) & (phase_gw < bins[i+1])
if np.any(mask):
mean_amp.append(np.mean(env_h1[mask]))
else:
mean_amp.append(0)
mean_amp = np.array(mean_amp)
if np.max(mean_amp) > 0:
mean_amp /= np.max(mean_amp)
plt.figure(figsize=(6,6))
plt.polar(bin_centers, mean_amp)
plt.title('Phase‑Amplitude Coupling: 200.2 Hz Envelope vs. GW Phase')
plt.grid(True)
plt.show()
mod_depth = (np.max(mean_amp) - np.min(mean_amp)) / np.mean(mean_amp) if np.mean(mean_amp) > 0 else 0
print(f"Modulation depth vs. GW phase: {mod_depth:.2f}\n")
# ---------------------------
# 8. FFT of envelope (look for orbital modulation)
# ---------------------------
# Use H1 envelope, avoid near‑merger edge
mask = t < -0.05
if np.sum(mask) > 100:
envelope_seg = env_h1[mask]
envelope_seg = envelope_seg - np.mean(envelope_seg)
fft_env = np.fft.rfft(envelope_seg)
freqs_env = np.fft.rfftfreq(len(envelope_seg), d=1/fs)
plt.figure(figsize=(10,4))
plt.semilogy(freqs_env, np.abs(fft_env))
plt.xlim(0, 200)
plt.xlabel('Modulation frequency (Hz)')
plt.ylabel('Power')
plt.title('Power Spectrum of 200.2 Hz Envelope (H1)')
f_min_orb = orbital_frequency(-0.45) # approx at middle of segment
f_max_orb = orbital_frequency(-0.05)
plt.axvspan(f_min_orb, f_max_orb, alpha=0.3, color='red', label=f'Expected f_orb range: {f_min_orb:.0f}–{f_max_orb:.0f} Hz')
plt.legend()
plt.grid(True)
plt.show()
peak_mod_idx = np.argmax(np.abs(fft_env)[1:]) + 1
peak_mod_freq = freqs_env[peak_mod_idx]
print(f"Dominant envelope modulation frequency: {peak_mod_freq:.1f} Hz")
mod_in_range = (f_min_orb <= peak_mod_freq <= f_max_orb)
print(f"Falls within expected orbital sweep? {'YES' if mod_in_range else 'NO'}\n")
else:
peak_mod_freq = np.nan
mod_in_range = False
print("Not enough envelope data for FFT.\n")
# ---------------------------
# 9. Final Summary
# ---------------------------
print("\n" + "="*60)
print("FINAL SUMMARY – FRCFD DEEP DIVE")
print("="*60)
print(f"1. Cross-correlation peak lag: {peak_lag_ms:.2f} ms (expected {DELAY_H1_L1*1000:.1f} ms) -> {'Match' if lag_match else 'No match'}")
print(f"2. Phase-amplitude modulation depth: {mod_depth:.2f} -> {'Significant (>0.2)' if mod_depth > 0.2 else 'Weak'}")
if not np.isnan(peak_mod_freq):
print(f"3. Envelope modulation frequency: {peak_mod_freq:.1f} Hz (expected orbital range {f_min_orb:.0f}–{f_max_orb:.0f} Hz) -> {'In range' if mod_in_range else 'Out of range'}")
else:
print("3. Envelope modulation frequency: Not computed (insufficient data)")
print("\nInterpretation:")
if lag_match and mod_depth > 0.2 and mod_in_range:
print("✅ All three indicators support the FRCFD regulator hypothesis.")
elif lag_match or mod_depth > 0.2 or mod_in_range:
print("⚠️ Partial evidence – some indicators agree, others do not.")
else:
print("❌ No clear evidence of FRCFD regulator in this dataset.")
print("="*60)
