import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt, find_peaks, welch

# --- Parameters ---
fs = 250                     # Sampling frequency (Hz)
duration_sec = 3600          # Duration (seconds)
n_samples = fs * duration_sec
time = np.arange(n_samples) / fs

HR_base = 75                 # Base heart rate (bpm)
# Increased modulation amplitudes for clearer PSD peaks
lf_freq, lf_amp = 0.1, 15    # LF frequency and amplitude (Hz, bpm)
hf_freq, hf_amp = 0.25, 7    # HF frequency and amplitude (Hz, bpm)

lowcut, highcut = 0.5, 40    # Bandpass filter cutoffs
filter_order = 4
window_size = 75             # Beats per HRV computation window (~1 minute)

# --- Helper functions ---

def bandpass_filter(signal, lowcut, highcut, fs, order=4):
    nyq = 0.5 * fs
    b, a = butter(order, [lowcut/nyq, highcut/nyq], btype='band')
    return filtfilt(b, a, signal)

def generate_hr(time):
    return HR_base + lf_amp * np.sin(2*np.pi*lf_freq*time) + hf_amp * np.sin(2*np.pi*hf_freq*time)

def generate_qrs_times(hr_inst, duration_sec, fs):
    rr_intervals = 60 / hr_inst
    qrs_times = []
    t, idx = 0, 0
    while t < duration_sec:
        interval = rr_intervals[min(idx, len(rr_intervals)-1)]
        t += interval
        qrs_times.append(t)
        idx += 1
    return np.array(qrs_times)

def synth_ecg(qrs_times, fs, length):
    ecg = np.zeros(length)
    waves = {
        'P':   {'amp': 0.2, 'offset': int(-0.2*fs), 'width': 0.02*fs},
        'QRS': {'amp': 1.0, 'offset': 0,             'width': 0.01*fs},
        'T':   {'amp': 0.3, 'offset': int(0.3*fs),  'width': 0.03*fs},
    }
    for tqrs in qrs_times:
        idx = int(tqrs * fs)
        for w in waves.values():
            c = idx + w['offset']
            start, end = c - int(w['width']*2), c + int(w['width']*2)
            if start < 0 or end >= length:
                continue
            xs = np.arange(start, end)
            ecg[start:end] += w['amp'] * np.exp(-0.5 * ((xs - c)/w['width'])**2)
    return ecg

def add_artifacts(ecg_signal, fs, fraction=0.0005, amplitude=3.0, duration_ms=10):
    ecg_art = ecg_signal.copy()
    n_spikes = int(len(ecg_signal) * fraction)
    spike_len = int(duration_ms * fs / 1000)
    positions = np.random.choice(len(ecg_signal) - spike_len, n_spikes, replace=False)
    for pos in positions:
        ecg_art[pos:pos+spike_len] += amplitude * np.hanning(spike_len)
    return ecg_art

def detect_qrs(ecg_signal, fs):
    filtered = bandpass_filter(ecg_signal, lowcut, highcut, fs, filter_order)
    norm_ecg = (filtered - np.mean(filtered)) / np.std(filtered)
    peaks, _ = find_peaks(norm_ecg, height=0.5, distance=int(0.4*fs))
    peak_times = peaks / fs
    return peak_times, filtered

def calc_hrv_time(rr_intervals, window_size):
    if len(rr_intervals) < window_size:
        return np.array([]), np.array([])
    rmssd, sdnn = [], []
    for i in range(0, len(rr_intervals)-window_size+1, window_size):
        seg = rr_intervals[i:i+window_size]
        diffs = np.diff(seg)
        rmssd.append(np.sqrt(np.mean(diffs**2)) if len(diffs)>0 else np.nan)
        sdnn.append(np.std(seg) if len(seg)>0 else np.nan)
    return np.array(rmssd), np.array(sdnn)

def calc_rr_psd(rr_intervals, fs_interp=8):
    import warnings
    if len(rr_intervals) < 2:
        return np.array([]), np.array([])
    time_cum = np.cumsum(rr_intervals)
    time_cum -= time_cum[0]  # Start at zero
    interp_time = np.arange(0, time_cum[-1], 1/fs_interp)
    if len(interp_time) < 2:
        return np.array([]), np.array([])
    interp_rr = np.interp(interp_time, time_cum, rr_intervals)
    interp_rr -= np.mean(interp_rr)  # Detrend RR interval signal
    nperseg_val = min(2048, len(interp_rr))
    if nperseg_val < 128:
        return np.array([]), np.array([])
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')
        f, psd = welch(interp_rr, fs=fs_interp, nperseg=nperseg_val)
    return f, psd

# --- Main pipeline execution ---

print("Step 1: Generate clean synthetic ECG...")
instant_hr = generate_hr(time)
qrs_times = generate_qrs_times(instant_hr, duration_sec, fs)
ecg_clean = synth_ecg(qrs_times, fs, n_samples)

plt.figure(figsize=(14,4))
plt.plot(time[:2500], ecg_clean[:2500])
plt.title("Step 1: Clean Synthetic ECG Signal (First 10 sec)")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude (mV)")
plt.grid(True)
plt.show()

print("Step 2: Add noise and artifacts...")
baseline_wander = 0.1 * np.sin(2*np.pi*0.05*time)
noise = 0.05 * np.random.randn(n_samples)
ecg_noisy = ecg_clean + baseline_wander + noise
ecg_noisy = add_artifacts(ecg_noisy, fs)

zoom_start = 1000
zoom_end = zoom_start + int(3*fs)
plt.figure(figsize=(14,4))
plt.plot(time[zoom_start:zoom_end], ecg_noisy[zoom_start:zoom_end])
plt.title("Step 2: Noisy Synthetic ECG with Artifacts (Zoomed View)")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude (mV)")
plt.grid(True)
plt.show()

print("Step 3: Filter ECG and detect QRS...")
detected_peaks, ecg_filtered = detect_qrs(ecg_noisy, fs)
print(f"Detected beats: {len(detected_peaks)}")

plt.figure(figsize=(14,4))
plt.plot(time[:2500], ecg_filtered[:2500], label='Filtered ECG')
plt.scatter(detected_peaks[detected_peaks < 10],
            np.interp(detected_peaks[detected_peaks < 10]*fs,
                      np.arange(len(ecg_filtered)), ecg_filtered),
            color='red', s=50, label='Detected R-peaks')
plt.title("Step 3: Filtered ECG and Detected R-peaks (First 10 sec)")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude (mV)")
plt.legend()
plt.grid(True)
plt.show()

print("Step 4: Calculate RR intervals and mean HR...")
rr_intervals = np.diff(detected_peaks)
print(f"Mean RR: {np.mean(rr_intervals):.4f}s, std RR: {np.std(rr_intervals):.6f}s")

plt.figure(figsize=(14,4))
plt.plot(detected_peaks[1:], rr_intervals*1000)
plt.title("Step 4: Beat-to-beat RR Intervals")
plt.xlabel("Time (s)")
plt.ylabel("RR Interval (ms)")
plt.grid(True)
plt.show()

mean_hr = 60 / np.mean(rr_intervals)
print(f"Mean Heart Rate: {mean_hr:.2f} bpm")

print("Step 5: Calculate time-domain HRV (RMSSD, SDNN)...")
rmssd, sdnn = calc_hrv_time(rr_intervals, window_size)

if len(rmssd) > 0:
    time_hrv = np.arange(len(rmssd)) * (window_size / HR_base)
    plt.figure(figsize=(14,4))
    plt.plot(time_hrv, rmssd, label='RMSSD')
    plt.plot(time_hrv, sdnn, label='SDNN')
    plt.title(f"Step 5: Time-domain HRV Metrics (Window size = {window_size})")
    plt.xlabel("Time (minutes)")
    plt.ylabel("Seconds")
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print("Not enough RR intervals for time-domain HRV.")

print("Step 6: Frequency-domain HRV analysis...")
f, psd = calc_rr_psd(rr_intervals)
if len(f) > 0:
    plt.figure(figsize=(14,4))
    plt.semilogy(f, psd, label='RR interval PSD')
    plt.axvline(lf_freq, color='orange', linestyle='--', label=f'LF ~{lf_freq} Hz')
    plt.axvline(hf_freq, color='green', linestyle='--', label=f'HF ~{hf_freq} Hz')
    plt.xlim(0, 0.5)
    plt.title("Step 6: Frequency domain HRV - PSD")
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Power Spectral Density (s^2/Hz)")
    plt.legend()
    plt.grid(True)
    plt.show()

    lf_band = (f >= 0.04) & (f <= 0.15)
    hf_band = (f >= 0.15) & (f <= 0.40)

    lf_power = np.trapezoid(psd[lf_band], f[lf_band]) if np.any(lf_band) else 0.0
    hf_power = np.trapezoid(psd[hf_band], f[hf_band]) if np.any(hf_band) else 0.0

    print(f"Low Frequency (LF) band power (0.04-0.15 Hz): {lf_power:.6f} s²/Hz")
    print(f"High Frequency (HF) band power (0.15-0.40 Hz): {hf_power:.6f} s²/Hz")
else:
    print("Not enough data for frequency domain HRV.")

print("\n*** SUMMARY ***")
print(f"Total detected beats: {len(detected_peaks)} (expected ~{len(qrs_times)})")
print(f"Mean Heart Rate: {mean_hr:.2f} bpm")
if len(f) > 0:
    print(f"LF power: {lf_power:.6f} s²/Hz, HF power: {hf_power:.6f} s²/Hz")
else:
    print("Frequency domain HRV metrics unavailable.")

print("ECG to HRV simulation complete.")


output

Low Frequency (LF) band power (0.04-0.15 Hz): 0.000000 s²/Hz
High Frequency (HF) band power (0.15-0.40 Hz): 0.000001 s²/Hz

*** SUMMARY ***
Total detected beats: 4482 (expected ~4482)
Mean Heart Rate: 74.70 bpm
LF power: 0.000000 s²/Hz, HF power: 0.000001 s²/Hz
ECG to HRV simulation complete.