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

# --- Parameters ---
# ECG sample frequency (Hz)
fs = 250
# 1 hour recording duration for synthetic data
duration_sec = 3600
# Total number of samples for the synthetic simulation
n_samples = fs * duration_sec
# Time array for the synthetic simulation
time = np.arange(n_samples) / fs

# Heart Rate (HR) control signals for Low Frequency (LF) and High Frequency (HF) modulations
HR_base = 75
lf_freq = 0.1  # Low-frequency modulation frequency (Hz)
hf_freq = 0.25 # High-frequency modulation frequency (Hz)
lf_amp = 5     # Amplitude of LF modulation
hf_amp = 2     # Amplitude of HF modulation

# ECG filter parameters
lowcut = 0.5
highcut = 40
filter_order = 4

# HRV window size for time-domain metrics (75 beats approximately 1-minute windows for an HR of 75 bpm)
window_size = 75

# ---------------------------------------------
# Helper functions

def bandpass_filter(signal, lowcut, highcut, fs, order=4):
    """
    Applies a Butterworth bandpass filter to the input signal.

    Args:
        signal (np.array): The input signal to filter.
        lowcut (float): The lower cutoff frequency of the bandpass filter.
        highcut (float): The upper cutoff frequency of the bandpass filter.
        fs (int): The sampling frequency of the signal.
        order (int, optional): The order of the filter. Defaults to 4.

    Returns:
        np.array: The filtered signal.
    """
    nyq = 0.5 * fs # Nyquist frequency
    # Calculate filter coefficients
    b, a = butter(order, [lowcut/nyq, highcut/nyq], btype='band')
    # Apply the filter forward and backward to avoid phase distortion
    return filtfilt(b, a, signal)

def generate_hr(time):
    """
    Generates instantaneous heart rate over time with LF and HF modulations.

    Args:
        time (np.array): Time array.

    Returns:
        np.array: Array of instantaneous heart rates.
    """
    # Base HR modulated by sine waves for LF and HF components
    hr = HR_base + lf_amp * np.sin(2*np.pi*lf_freq*time) + hf_amp * np.sin(2*np.pi*hf_freq*time)
    return hr

def generate_qrs_times(hr_inst, duration_sec, fs):
    """
    From instantaneous HR, generates beat times as a cumulative sum of RR intervals.

    Args:
        hr_inst (np.array): Instantaneous heart rate array.
        duration_sec (int): Total duration of the recording in seconds.
        fs (int): Sampling frequency.

    Returns:
        np.array: Array of QRS peak times in seconds.
    """
    # Convert instantaneous HR to RR intervals (seconds)
    rr_intervals = 60 / hr_inst
    qrs_times = []
    t = 0
    idx = 0
    # Iterate to generate cumulative QRS times
    while t < duration_sec:
        # Get the current RR interval, using the last one if out of bounds
        interval = rr_intervals[idx] if idx < len(rr_intervals) else rr_intervals[-1]
        t += interval # Accumulate time for the next beat
        qrs_times.append(t)
        idx += 1
    return np.array(qrs_times)

def synth_ecg(qrs_times, fs, length):
    """
    Generates a synthetic ECG signal with P, QRS, and T waves.

    Args:
        qrs_times (np.array): Times of QRS peaks in seconds.
        fs (int): Sampling frequency.
        length (int): Total number of samples for the ECG signal.

    Returns:
        np.array: Synthetic ECG signal.
    """
    ecg = np.zeros(length)
    # Define typical durations/amplitudes for P, QRS, T waves relative to sampling frequency
    # QRS complex (main peak)
    qrs_peak_amplitude = 1.0 # mV
    qrs_width_factor = 0.01 * fs # controls sharpness of QRS

    # P wave (atrial depolarization)
    p_wave_amplitude = 0.2 # mV
    p_wave_offset_samples = int(-0.2 * fs) # P wave occurs ~200ms before R-peak
    p_wave_width_factor = 0.02 * fs

    # T wave (ventricular repolarization)
    t_wave_amplitude = 0.3 # mV
    t_wave_offset_samples = int(0.3 * fs) # T wave occurs ~300ms after R-peak
    t_wave_width_factor = 0.03 * fs

    for tqrs in qrs_times:
        # Convert QRS time to sample index
        idx = int(tqrs * fs)

        # Generate QRS complex (Gaussian-shaped)
        qrs_indices = np.arange(idx - int(0.05 * fs), idx + int(0.05 * fs))
        valid_qrs = (qrs_indices >= 0) & (qrs_indices < length)
        ecg[qrs_indices[valid_qrs]] += qrs_peak_amplitude * np.exp(-0.5 * ((qrs_indices[valid_qrs] - idx)/qrs_width_factor)**2)

        # Generate P wave (Gaussian-shaped)
        p_idx = idx + p_wave_offset_samples
        p_indices = np.arange(p_idx - int(0.05 * fs), p_idx + int(0.05 * fs))
        valid_p = (p_indices >= 0) & (p_indices < length)
        ecg[p_indices[valid_p]] += p_wave_amplitude * np.exp(-0.5 * ((p_indices[valid_p] - p_idx)/p_wave_width_factor)**2)

        # Generate T wave (Gaussian-shaped)
        t_idx = idx + t_wave_offset_samples
        t_indices = np.arange(t_idx - int(0.08 * fs), t_idx + int(0.08 * fs))
        valid_t = (t_indices >= 0) & (t_indices < length)
        ecg[t_indices[valid_t]] += t_wave_amplitude * np.exp(-0.5 * ((t_indices[valid_t] - t_idx)/t_wave_width_factor)**2)

    return ecg

def add_artifacts(ecg_signal, fs, fraction=0.001, amplitude=3.0, duration_ms=10):
    """
    Adds spike artifacts to the ECG signal, simulating electrode pops or muscle noise.

    Args:
        ecg_signal (np.array): The clean ECG signal.
        fs (int): Sampling frequency.
        fraction (float, optional): Fraction of the signal length where spikes will occur. Defaults to 0.001.
        amplitude (float, optional): Amplitude of the spikes. Defaults to 3.0.
        duration_ms (int, optional): Duration of each spike in milliseconds. Defaults to 10.

    Returns:
        np.array: ECG signal with added artifacts.
    """
    ecg_art = ecg_signal.copy()
    # Number of spikes to add
    n_spikes = int(len(ecg_signal) * fraction)
    # Duration of a spike in samples
    spike_duration_samples = int(duration_ms * fs / 1000)
    # Randomly choose positions for the spikes
    spike_positions = np.random.choice(len(ecg_signal) - spike_duration_samples, n_spikes, replace=False)
    for pos in spike_positions:
        # Add a Hanning window shaped spike
        ecg_art[pos:pos + spike_duration_samples] += amplitude * np.hanning(spike_duration_samples)
    return ecg_art

def detect_qrs(ecg_signal, fs):
    """
    Detects QRS peaks using simple thresholding on a filtered signal.

    Args:
        ecg_signal (np.array): The input ECG signal (potentially noisy).
        fs (int): Sampling frequency.

    Returns:
        tuple: A tuple containing:
            - np.array: Times of detected QRS peaks in seconds.
            - np.array: The filtered ECG signal.
    """
    # Filter the ECG signal to enhance QRS complexes
    filtered = bandpass_filter(ecg_signal, lowcut, highcut, fs, filter_order)
    # Normalize the filtered signal
    norm_ecg = (filtered - np.mean(filtered)) / np.std(filtered)
    # Find peaks using `find_peaks` with a height threshold and minimum distance
    peaks, _ = find_peaks(norm_ecg, height=0.5, distance=int(0.4*fs)) # 0.4s refractory period
    peak_times = peaks / fs # Convert peak indices to times
    return peak_times, filtered

def calc_hrv_time(rrs, window_size):
    """
    Calculates time-domain HRV metrics (RMSSD and SDNN) over sliding windows.

    Args:
        rrs (np.array): Array of RR intervals.
        window_size (int): The number of beats in each window.

    Returns:
        tuple: A tuple containing:
            - np.array: RMSSD values for each window.
            - np.array: SDNN values for each window.
    """
    rmssd = []
    sdnn = []
    # Iterate through the RR intervals with a sliding window
    # Ensure there are enough beats for at least one window
    if len(rrs) < window_size:
        return np.array([]), np.array([])

    for i in range(0, len(rrs) - window_size + 1, window_size):
        seg = rrs[i:i+window_size] # Current window segment
        # Calculate RMSSD (Root Mean Square of Successive Differences)
        diffs = np.diff(seg)
        rmssd.append(np.sqrt(np.mean(diffs**2)) if len(diffs) > 0 else np.nan)
        # Calculate SDNN (Standard Deviation of NN intervals)
        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=4):
    """
    Calculates the Power Spectral Density (PSD) of RR intervals for frequency-domain HRV.

    Args:
        rr_intervals (np.array): Array of RR intervals.
        fs_interp (int, optional): Interpolation frequency for resampling RR intervals. Defaults to 4 Hz.

    Returns:
        tuple: A tuple containing:
            - np.array: Frequencies corresponding to the PSD.
            - np.array: Power Spectral Density values.
    """
    if len(rr_intervals) < 2:
        return np.array([]), np.array([]) # Not enough data for PSD

    # Calculate cumulative time from RR intervals
    time_cum = np.cumsum(rr_intervals)
    # Ensure interp_time spans a valid range
    if time_cum[-1] - time_cum[0] <= 0:
        return np.array([]), np.array([])

    # Create an interpolated time axis
    interp_time = np.arange(time_cum[0], time_cum[-1], 1/fs_interp)
    if len(interp_time) < 2: # Ensure at least 2 points for interpolation
        return np.array([]), np.array([])

    # Interpolate RR intervals onto the new time axis
    interp_rr = np.interp(interp_time, time_cum, rr_intervals)
    interp_rr -= np.mean(interp_rr) # Detrend the interpolated RR series

    # Calculate PSD using Welch's method
    # nperseg should be <= len(interp_rr)
    nperseg_val = min(1024, len(interp_rr))
    if nperseg_val < 2: # Need at least 2 points for welch
        return np.array([]), np.array([])

    f, psd = welch(interp_rr, fs=fs_interp, nperseg=nperseg_val)
    return f, psd


# ---------------------------------------------
# STEP 1: Generate Raw Clean Synthetic ECG signal (with P, QRS, T waves)

print("Step 1: Generating raw clean synthetic ECG signal with P, QRS, and T waves...")
# Generate instantaneous heart rate profile
instant_hr = generate_hr(time)
# Generate QRS peak times based on the instantaneous HR
qrs_times_synth = generate_qrs_times(instant_hr, duration_sec, fs)
# Synthesize the clean ECG signal with more realistic morphology
ecg_clean_synth = synth_ecg(qrs_times_synth, fs, n_samples)

# Plot the first 10 seconds of the clean synthetic ECG signal
plt.figure(figsize=(14,4))
plt.plot(time[:2500], ecg_clean_synth[:2500])
plt.title("Step 1: Clean Synthetic ECG Signal (First 10 seconds) - Ideal Heartbeats")
plt.xlabel("Time (seconds)")
plt.ylabel("Amplitude (mV)")
plt.grid(True)
plt.show()

# ---------------------------------------------
# STEP 2: Add noise (baseline wander, Gaussian noise, and spike artifacts) to synthetic ECG

print("Step 2: Adding realistic noise and spike artifacts to synthetic ECG...")
# Simulate baseline wander (slow sinusoidal drift)
baseline_wander = 0.1 * np.sin(2*np.pi*0.05*time)
# Simulate high-frequency Gaussian noise
noise = 0.05 * np.random.randn(n_samples)
# Combine clean ECG with baseline wander and Gaussian noise
ecg_with_noise_synth = ecg_clean_synth + baseline_wander + noise
# Add transient spike artifacts to the noisy ECG signal
ecg_realistic_noise_synth = add_artifacts(ecg_with_noise_synth, fs, fraction=0.0005, amplitude=3.0, duration_ms=10)

# Plot a zoomed-in section to show the combined noise and artifacts clearly
zoom_start_synth = 1000 # Start index for zooming (approx 4 seconds)
zoom_end_synth = zoom_start_synth + int(3*fs) # End index for zooming (3 seconds window)

plt.figure(figsize=(14,4))
plt.plot(time[zoom_start_synth:zoom_end_synth], ecg_realistic_noise_synth[zoom_start_synth:zoom_end_synth], label="Noisy ECG Signal")
plt.title("Step 2: Synthetic ECG Signal with Added Noise and Artifacts (Zoomed View)")
plt.xlabel("Time (seconds)")
plt.ylabel("Amplitude (mV)")
plt.grid(True)
plt.show()

# ---------------------------------------------
# STEP 3: Signal processing and QRS detection on Synthetic ECG

print("Step 3: Filtering synthetic ECG and detecting QRS peaks...")
detected_qrs_times_synth, ecg_filtered_synth = detect_qrs(ecg_realistic_noise_synth, fs)
print(f"Synthetic ECG - Number of beats detected: {len(detected_qrs_times_synth)} beats")

# Plot the first 10 seconds of the filtered synthetic ECG with detected QRS peaks
plt.figure(figsize=(14,4))
plt.plot(time[:2500], ecg_filtered_synth[:2500], label="Filtered ECG Signal")
plt.scatter(detected_qrs_times_synth[detected_qrs_times_synth < 10],
            np.interp(detected_qrs_times_synth[detected_qrs_times_synth < 10]*fs, np.arange(len(ecg_filtered_synth)), ecg_filtered_synth[0:]),
            color='red', marker='o', s=50, label="Detected R-peaks")
plt.title("Step 3.1: Filtered Synthetic ECG and Detected R-peaks (First 10 seconds)")
plt.xlabel("Time (seconds)")
plt.ylabel("Amplitude (mV)")
plt.legend()
plt.grid(True)
plt.show()

# ---------------------------------------------
# STEP 4: RR interval calculation & plotting (from Synthetic ECG for the full duration)

print("Step 4: Calculate RR intervals from detected beats (Synthetic ECG, full 1 hour)...")
rr_intervals_synth = np.diff(detected_qrs_times_synth)
mean_hr_synth = 60 / np.mean(rr_intervals_synth) if len(rr_intervals_synth) > 0 else np.nan
print(f"Synthetic ECG - Mean heart rate: {mean_hr_synth:.2f} bpm")

# Plot RR intervals over 1 hour (Synthetic)
plt.figure(figsize=(14,4))
plt.plot(detected_qrs_times_synth[1:], rr_intervals_synth*1000, label="RR Intervals") # Convert RR to milliseconds
plt.title("Step 4.1: Synthetic ECG - Beat-to-Beat RR Intervals over 1 hour")
plt.xlabel("Time (seconds)")
plt.ylabel("RR Interval (milliseconds)")
plt.grid(True)
plt.show()

# ---------------------------------------------
# STEP 5: Calculate and plot time-domain HRV metrics (Synthetic)

print("Step 5: Calculate time-domain HRV metrics...")
rmssd_synth, sdnn_synth = calc_hrv_time(rr_intervals_synth, window_size)

if len(rmssd_synth) > 0:
    window_duration_min = window_size / HR_base # Approximate window duration in minutes
    x_time_synth_hrv = np.arange(len(rmssd_synth)) * window_duration_min

    plt.figure(figsize=(14,4))
    plt.plot(x_time_synth_hrv, rmssd_synth, label="RMSSD (Root Mean Square of Successive Differences)")
    plt.plot(x_time_synth_hrv, sdnn_synth, label="SDNN (Standard Deviation of NN Intervals)")
    plt.title(f"Step 5.1: Synthetic ECG - Time-domain HRV Metrics (Window size = {window_size} beats)")
    plt.xlabel("Time (minutes)")
    plt.ylabel("HRV Metric Value (seconds)")
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print("Not enough synthetic RR intervals to calculate time-domain HRV metrics.")


# ---------------------------------------------
# STEP 6: Calculate and plot frequency-domain HRV (Synthetic)

print("Step 6: Calculate frequency domain PSD...")
f_synth, psd_synth = calc_rr_psd(rr_intervals_synth)

if len(f_synth) > 0:
    plt.figure(figsize=(14,4))
    plt.semilogy(f_synth, psd_synth, label="Power Spectral Density of RR Intervals")
    plt.axvline(lf_freq, color='orange', linestyle='--', label=f"LF Band (~{lf_freq} Hz)")
    plt.axvline(hf_freq, color='green', linestyle='--', label=f"HF Band (~{hf_freq} Hz)")
    plt.title("Step 6.1: Synthetic ECG - Frequency Domain HRV - Power Spectral Density")
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Power Spectral Density (s²/Hz)")
    plt.xlim(0, 0.5) # Limit x-axis to relevant HRV frequency range
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print("Not enough synthetic RR intervals to calculate frequency-domain HRV metrics.")

# ---------------------------------------------
# Final summary

print("\n*** SUMMARY ***")
print("\n--- Synthetic ECG Pipeline Summary ---")
print(f"Total beats detected: {len(detected_qrs_times_synth)} (expected ~{len(qrs_times_synth)})")
print(f"Mean Heart Rate from detected beats: {mean_hr_synth:.2f} bpm")
if len(f_synth) > 0:
    # Changed np.trapz to np.trapezoid as per deprecation warning
    lf_band_power_synth = np.trapezoid(psd_synth[(f_synth >= 0.04) & (f_synth <= 0.15)], f_synth[(f_synth >= 0.04) & (f_synth <= 0.15)])
    hf_band_power_synth = np.trapezoid(psd_synth[(f_synth >= 0.15) & (f_synth <= 0.40)], f_synth[(f_synth >= 0.15) & (f_synth <= 0.40)])
    print(f"Low Frequency (LF) band power (0.04-0.15 Hz): {lf_band_power_synth:.6f} s^2/Hz")
    print(f"High Frequency (HF) band power (0.15-0.40 Hz): {hf_band_power_synth:.6f} s^2/Hz")
else:
    print("Frequency domain HRV metrics not calculated for synthetic data due to insufficient data.")

print("\nThe comprehensive ECG-to-HRV analysis pipeline simulation is complete.")
print("The output simulates a wearable ECG patch's data generation and HRV analysis.")






this is final.. i think its more then enough for my task right?