Wearable ECG Patch - simulation -Heart rate variability (HRV) over time over 1h


The goal is to setup a signal processing simulation


Random values are ok for noise, as noise is filtered.


The goal is to validate your simulation at a certain point, thus you

must check if the result is correct. Using random values at the first

place, does not help. But you may add them later.

Regards,
M. Wiehl

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

np.random.seed(42)  # For reproducible results

def generate_rr_synthetic(num_beats=4500, base_hr=75, lf_amp=0.05, hf_amp=0.02,
                          lf_freq=0.1, hf_freq=0.25, fs_interp=4, noise_std=0.0085):
    """
    Generate synthetic RR intervals with LF and HF sinusoidal modulations plus Gaussian noise.
    """
    base_rr = 60 / base_hr
    t = np.arange(num_beats) / fs_interp
    lf = lf_amp * np.sin(2 * np.pi * lf_freq * t)
    hf = hf_amp * np.sin(2 * np.pi * hf_freq * t)
    noise = np.random.normal(0, noise_std, num_beats)
    rr = base_rr + lf + hf + noise
    return rr

def add_spikes(rr, spike_fraction=0.02, spike_magnitude=0.15):
    """
    Add artificial spikes to RR intervals to simulate ectopic beats or artifacts.
    """
    rr_spiked = rr.copy()
    num_spikes = int(len(rr) * spike_fraction)
    spike_indices = np.random.choice(len(rr), num_spikes, replace=False)
    random_signs = np.random.choice([-1, 1], num_spikes)
    spikes = random_signs * spike_magnitude * np.random.rand(num_spikes)
    rr_spiked[spike_indices] += spikes
    return rr_spiked

def calc_hrv_time_domain(rr, window_size=75):
    """
    Calculate RMSSD and SDNN in non-overlapping windows.
    """
    rmssd_vals = []
    sdnn_vals = []
    for start in range(0, len(rr) - window_size + 1, window_size):
        segment = rr[start:start+window_size]
        diffs = np.diff(segment)
        rmssd_vals.append(np.sqrt(np.mean(diffs**2)))
        sdnn_vals.append(np.std(segment))
    return np.array(rmssd_vals), np.array(sdnn_vals)

def rr_to_psd(rr, fs=4):
    """
    Calculate PSD of RR intervals interpolated at uniform sampling rate fs.
    Returns frequencies, PSD, LF power, and HF power.
    """
    time_cum = np.cumsum(rr)
    interp_time = np.arange(time_cum[0], time_cum[-1], 1/fs)
    interp_rr = np.interp(interp_time, time_cum, rr)
    interp_rr -= np.mean(interp_rr)
    f, psd = welch(interp_rr, fs=fs, nperseg=1024)
    lf_band = (f >= 0.04) & (f <= 0.15)
    hf_band = (f >= 0.15) & (f <= 0.4)
    lf_power = np.trapezoid(psd[lf_band], f[lf_band])
    hf_power = np.trapezoid(psd[hf_band], f[hf_band])
    return f, psd, lf_power, hf_power

def load_rr_from_file(filename, max_beats=4500):
    """
    Load RR intervals in ms from file and convert to seconds.
    """
    with open(filename, "r") as f:
        rr_ms = [line.strip() for line in f if line.strip().isdigit()]
    rr_ms = [int(x) for x in rr_ms]
    rr_sec = np.array(rr_ms) / 1000.0
    return rr_sec[:max_beats]

def plot_and_save():
    # Synthetic RR with spikes
    rr_syn = generate_rr_synthetic()
    rr_syn = add_spikes(rr_syn)

    rmssd_syn, sdnn_syn = calc_hrv_time_domain(rr_syn)
    f_syn, psd_syn, lf_syn, hf_syn = rr_to_psd(rr_syn)

    print(f"Synthetic LF power: {lf_syn:.5f} s²/Hz")
    print(f"Synthetic HF power: {hf_syn:.5f} s²/Hz")

    plt.figure(figsize=(16, 12))

    plt.subplot(3,1,1)
    plt.plot(rr_syn)
    plt.title("Synthetic RR Intervals (with spikes)")
    plt.xlabel("Beat Number")
    plt.ylabel("RR Interval (s)")
    plt.grid()

    plt.subplot(3,1,2)
    plt.plot(rmssd_syn, label="RMSSD")
    plt.plot(sdnn_syn, label="SDNN")
    plt.title("Time-domain HRV Metrics (Synthetic)")
    plt.xlabel("Window Number")
    plt.ylabel("Seconds")
    plt.legend()
    plt.grid()

    plt.subplot(3,1,3)
    plt.semilogy(f_syn, psd_syn)
    plt.axvline(0.1, color='orange', linestyle='--', label="LF Band ~0.1 Hz")
    plt.axvline(0.25, color='green', linestyle='--', label="HF Band ~0.25 Hz")
    plt.title("Frequency-domain HRV PSD (Synthetic)")
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Power spectral density (s²/Hz)")
    plt.xlim(0, 0.5)
    plt.legend()
    plt.grid()

    plt.tight_layout()
    plt.savefig("synthetic_hrv_analysis.png")
    plt.show()

    # Load and analyze real RR data
    try:
        rr_real = load_rr_from_file("000.txt")
        rmssd_real, sdnn_real = calc_hrv_time_domain(rr_real)
        f_real, psd_real, lf_real, hf_real = rr_to_psd(rr_real)

        print(f"Real LF power: {lf_real:.5f} s²/Hz")
        print(f"Real HF power: {hf_real:.5f} s²/Hz")

        plt.figure(figsize=(14,6))
        plt.plot(rmssd_syn[:len(rmssd_real)], label="Synthetic RMSSD")
        plt.plot(rmssd_real, label="Real RMSSD")
        plt.title("RMSSD: Synthetic vs Real")
        plt.xlabel("Window Number")
        plt.ylabel("Seconds")
        plt.legend()
        plt.grid()
        plt.savefig("rmssd_comparison.png")
        plt.show()

        plt.figure(figsize=(14,6))
        plt.plot(sdnn_syn[:len(sdnn_real)], label="Synthetic SDNN")
        plt.plot(sdnn_real, label="Real SDNN")
        plt.title("SDNN: Synthetic vs Real")
        plt.xlabel("Window Number")
        plt.ylabel("Seconds")
        plt.legend()
        plt.grid()
        plt.savefig("sdnn_comparison.png")
        plt.show()

        plt.figure(figsize=(14,6))
        plt.semilogy(f_syn, psd_syn, label="Synthetic PSD")
        plt.semilogy(f_real, psd_real, label="Real PSD")
        plt.title("PSD: Synthetic vs Real")
        plt.xlabel("Frequency (Hz)")
        plt.ylabel("Power spectral density (s²/Hz)")
        plt.xlim(0, 0.5)
        plt.legend()
        plt.grid()
        plt.savefig("psd_comparison.png")
        plt.show()

    except FileNotFoundError:
        print("Real data file '000.txt' not found. Skipping real data analysis.")
        
def print_statistical_summary(label, rmssd, sdnn, lf_power, hf_power):
    print(f"Statistical Summary for {label}:")
    print(f"{'Metric':<10} {'Mean':>10} {'Std Dev':>10}")
    print(f"{'-'*32}")
    print(f"{'RMSSD (s)':<10} {np.mean(rmssd):10.4f} {np.std(rmssd):10.4f}")
    print(f"{'SDNN (s)':<10} {np.mean(sdnn):10.4f} {np.std(sdnn):10.4f}")
    print(f"{'LF Power':<10} {lf_power:10.6f} {'N/A':>10}")
    print(f"{'HF Power':<10} {hf_power:10.6f} {'N/A':>10}")
    print()

def simulate_healthy_vs_stressed():
    # Healthy: lower HR, higher variability
    rr_healthy = generate_rr_synthetic(base_hr=70, lf_amp=0.05, hf_amp=0.02)
    rr_healthy = add_spikes(rr_healthy, spike_fraction=0.02, spike_magnitude=0.15)

    # Stressed: higher HR, lower variability
    rr_stressed = generate_rr_synthetic(base_hr=90, lf_amp=0.02, hf_amp=0.01)
    rr_stressed = add_spikes(rr_stressed, spike_fraction=0.02, spike_magnitude=0.15)

    rmssd_healthy, sdnn_healthy = calc_hrv_time_domain(rr_healthy)
    rmssd_stressed, sdnn_stressed = calc_hrv_time_domain(rr_stressed)

    plt.figure(figsize=(12,5))

    plt.subplot(1,2,1)
    plt.plot(rmssd_healthy, label="Healthy", color='green')
    plt.plot(rmssd_stressed, label="Stressed", color='red', linestyle='--')
    plt.title("RMSSD: Healthy vs Stressed")
    plt.xlabel("Window Number")
    plt.ylabel("Seconds")
    plt.legend()
    plt.grid()

    plt.subplot(1,2,2)
    plt.plot(sdnn_healthy, label="Healthy", color='green')
    plt.plot(sdnn_stressed, label="Stressed", color='red', linestyle='--')
    plt.title("SDNN: Healthy vs Stressed")
    plt.xlabel("Window Number")
    plt.ylabel("Seconds")
    plt.legend()
    plt.grid()

    plt.tight_layout()
    plt.savefig("healthy_vs_stressed_hrv.png")
    plt.show()

if __name__ == "__main__":
    plot_and_save()  # Run synthetic + real data analysis
    simulate_healthy_vs_stressed()  # Run healthy vs stressed comparison






is it according to intruction?

as i have to simulate ECG patch like what data the ECG patch generate while working RAW one and then that singnal how it process and give usefull data and show the hrv over time over 1 hour.. 
whole process i have to simulate.. as real ECG patch might works..