Tools: Coding the Heart: Using SciPy and HRV to Predict Your Burnout Level

The Architecture of Heart Rate Variability (HRV) Analysis ## Prerequisites ## Step 1: Denoising the Signal with a Kalman Filter ## Step 2: Extracting R-Peaks and NN-Intervals ## Step 3: Frequency Domain Analysis (The LF/HF Ratio) ## The "Official" Way to Scale ## Conclusion: Listen to Your Data As developers, we are experts at debugging complex microservices, but we often ignore the most important system we'll ever maintain: our own body. Have you ever felt "fried" after a long sprint? That’s not just mental fatigue; it’s your Autonomic Nervous System (ANS) signaling for help. In this tutorial, we will dive deep into HRV signal processing and frequency domain analysis. We'll learn how to transform raw wearable data into actionable insights using the Python scientific stack. By calculating the LF/HF ratio, we can quantitatively assess stress levels and stop burnout before the "system crash" happens. We’ll be utilizing powerful tools like SciPy, NumPy, and BioSPPy to handle the heavy lifting of signal processing. If you’re interested in more production-ready patterns for health-tech and signal processing, I highly recommend checking out the advanced case studies at WellAlly Tech Blog, which served as a major inspiration for this technical deep dive. To go from a noisy voltage signal (ECG/PPG) to a burnout index, we follow a specific signal-processing pipeline. Here is how the data flows from your wrist to your dashboard: Ensure you have your environment ready. We will use BioSPPy for biosignal processing and SciPy for the mathematical transformations. Wearable sensors are notoriously noisy due to "motion artifacts" (moving your arm while typing). A standard moving average isn't enough; we need a Kalman Filter to estimate the true state of the signal. The "R-peak" is the highest point in the heartbeat. The time between these peaks (NN-intervals) is what we actually measure for HRV. We'll use the BioSPPy library, which provides a robust implementation of the Hamilton segmenter. This is where the magic happens. We decompose the heart rhythm into frequency components: Burnout Metric: A high LF/HF ratio indicates sympathetic dominance—meaning you are stressed, caffeinated, or overworked. While the script above works for local analysis, building a production-grade health monitoring system requires handling edge cases like ectopic beats, signal dropouts, and real-time streaming. If you want to see how these algorithms are integrated into large-scale cloud architectures or how to optimize signal processing for mobile chips, check out the deep dives at wellally.tech/blog. They cover the intersection of AI, signal processing, and high-performance computing in the healthcare space—essential reading for any dev moving into the Bio-IT sector. By monitoring your LF/HF ratio, you can see burnout coming before you start making "fatigue-driven" commits. A ratio consistently above 2.0 during rest is a clear sign that you need to step away from the keyboard, hydrate, and maybe—just maybe—close those 47 open Chrome tabs. What's your strategy for tracking health metrics? Have you tried analyzing your Oura or Apple Watch data? Drop a comment below! 🚀 Templates let you quickly answer FAQs or store snippets for re-use. Are you sure you want to hide this comment? It will become hidden in your post, but will still be visible via the comment's permalink. Hide child comments as well For further actions, you may consider blocking this person and/or reporting abuse COMMAND_BLOCK: graph TD A[Raw ECG/PPG Signal] --> B[Denoising: Kalman Filter] B --> C[R-Peak Detection: BioSPPy] C --> D[NN-Interval Calculation] D --> E[Resampling & Interpolation] E --> F[Frequency Domain Analysis: FFT/Welch] F --> G[LF/HF Ratio calculation] G --> H[Burnout Index Assessment] Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: graph TD A[Raw ECG/PPG Signal] --> B[Denoising: Kalman Filter] B --> C[R-Peak Detection: BioSPPy] C --> D[NN-Interval Calculation] D --> E[Resampling & Interpolation] E --> F[Frequency Domain Analysis: FFT/Welch] F --> G[LF/HF Ratio calculation] G --> H[Burnout Index Assessment] COMMAND_BLOCK: graph TD A[Raw ECG/PPG Signal] --> B[Denoising: Kalman Filter] B --> C[R-Peak Detection: BioSPPy] C --> D[NN-Interval Calculation] D --> E[Resampling & Interpolation] E --> F[Frequency Domain Analysis: FFT/Welch] F --> G[LF/HF Ratio calculation] G --> H[Burnout Index Assessment] COMMAND_BLOCK: pip install numpy scipy matplotlib biosppy Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: pip install numpy scipy matplotlib biosppy COMMAND_BLOCK: pip install numpy scipy matplotlib biosppy COMMAND_BLOCK: import numpy as np from scipy.signal import standard_scaler def kalman_filter(signal): """ A simple 1D Kalman Filter to smooth raw ECG data. """ n_iter = len(signal) sz = (n_iter,) z = signal # Observations Q = 1e-5 # Process variance R = 0.1**2 # Measurement variance xhat = np.zeros(sz) # Posteriori estimate P = np.zeros(sz) # Posteriori error covariance xhatminus = np.zeros(sz) # Priori estimate Pminus = np.zeros(sz) # Priori error covariance K = np.zeros(sz) # Kalman gain xhat[0] = z[0] P[0] = 1.0 for k in range(1, n_iter): # Time Update xhatminus[k] = xhat[k-1] Pminus[k] = P[k-1] + Q # Measurement Update K[k] = Pminus[k] / (Pminus[k] + R) xhat[k] = xhatminus[k] + K[k] * (z[k] - xhatminus[k]) P[k] = (1 - K[k]) * Pminus[k] return xhat Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: import numpy as np from scipy.signal import standard_scaler def kalman_filter(signal): """ A simple 1D Kalman Filter to smooth raw ECG data. """ n_iter = len(signal) sz = (n_iter,) z = signal # Observations Q = 1e-5 # Process variance R = 0.1**2 # Measurement variance xhat = np.zeros(sz) # Posteriori estimate P = np.zeros(sz) # Posteriori error covariance xhatminus = np.zeros(sz) # Priori estimate Pminus = np.zeros(sz) # Priori error covariance K = np.zeros(sz) # Kalman gain xhat[0] = z[0] P[0] = 1.0 for k in range(1, n_iter): # Time Update xhatminus[k] = xhat[k-1] Pminus[k] = P[k-1] + Q # Measurement Update K[k] = Pminus[k] / (Pminus[k] + R) xhat[k] = xhatminus[k] + K[k] * (z[k] - xhatminus[k]) P[k] = (1 - K[k]) * Pminus[k] return xhat COMMAND_BLOCK: import numpy as np from scipy.signal import standard_scaler def kalman_filter(signal): """ A simple 1D Kalman Filter to smooth raw ECG data. """ n_iter = len(signal) sz = (n_iter,) z = signal # Observations Q = 1e-5 # Process variance R = 0.1**2 # Measurement variance xhat = np.zeros(sz) # Posteriori estimate P = np.zeros(sz) # Posteriori error covariance xhatminus = np.zeros(sz) # Priori estimate Pminus = np.zeros(sz) # Priori error covariance K = np.zeros(sz) # Kalman gain xhat[0] = z[0] P[0] = 1.0 for k in range(1, n_iter): # Time Update xhatminus[k] = xhat[k-1] Pminus[k] = P[k-1] + Q # Measurement Update K[k] = Pminus[k] / (Pminus[k] + R) xhat[k] = xhatminus[k] + K[k] * (z[k] - xhatminus[k]) P[k] = (1 - K[k]) * Pminus[k] return xhat COMMAND_BLOCK: from biosppy.signals import ecg def get_nn_intervals(raw_signal, sampling_rate=1000): # Process the ECG signal # This automatically filters and finds R-peaks out = ecg.ecg(signal=raw_signal, sampling_rate=sampling_rate, show=False) # R-peak indices rpeaks = out['rpeaks'] # Calculate intervals in milliseconds nn_intervals = np.diff(rpeaks) * (1000.0 / sampling_rate) return nn_intervals Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: from biosppy.signals import ecg def get_nn_intervals(raw_signal, sampling_rate=1000): # Process the ECG signal # This automatically filters and finds R-peaks out = ecg.ecg(signal=raw_signal, sampling_rate=sampling_rate, show=False) # R-peak indices rpeaks = out['rpeaks'] # Calculate intervals in milliseconds nn_intervals = np.diff(rpeaks) * (1000.0 / sampling_rate) return nn_intervals COMMAND_BLOCK: from biosppy.signals import ecg def get_nn_intervals(raw_signal, sampling_rate=1000): # Process the ECG signal # This automatically filters and finds R-peaks out = ecg.ecg(signal=raw_signal, sampling_rate=sampling_rate, show=False) # R-peak indices rpeaks = out['rpeaks'] # Calculate intervals in milliseconds nn_intervals = np.diff(rpeaks) * (1000.0 / sampling_rate) return nn_intervals COMMAND_BLOCK: from scipy.interpolate import interp1d from scipy.signal import welch def analyze_frequency_domain(nn_intervals): # 1. Resample: NN-intervals are irregularly spaced. We need 4Hz interpolation for FFT. time = np.cumsum(nn_intervals) / 1000.0 # seconds f_interp = interp1d(time, nn_intervals, kind='cubic') new_time = np.arange(time[0], time[-1], 0.25) # 4Hz nn_resampled = f_interp(new_time) # 2. Power Spectral Density using Welch's method freqs, psd = welch(nn_resampled, fs=4.0, nperseg=256) # 3. Define bands lf_band = (freqs >= 0.04) & (freqs <= 0.15) hf_band = (freqs >= 0.15) & (freqs <= 0.4) # 4. Calculate Area Under Curve (AUC) lf_power = np.trapz(psd[lf_band], freqs[lf_band]) hf_power = np.trapz(psd[hf_band], freqs[hf_band]) return lf_power / hf_power # Usage # ratio = analyze_frequency_domain(my_intervals) # print(f"Burnout Index (LF/HF): {ratio:.2f}") Enter fullscreen mode Exit fullscreen mode COMMAND_BLOCK: from scipy.interpolate import interp1d from scipy.signal import welch def analyze_frequency_domain(nn_intervals): # 1. Resample: NN-intervals are irregularly spaced. We need 4Hz interpolation for FFT. time = np.cumsum(nn_intervals) / 1000.0 # seconds f_interp = interp1d(time, nn_intervals, kind='cubic') new_time = np.arange(time[0], time[-1], 0.25) # 4Hz nn_resampled = f_interp(new_time) # 2. Power Spectral Density using Welch's method freqs, psd = welch(nn_resampled, fs=4.0, nperseg=256) # 3. Define bands lf_band = (freqs >= 0.04) & (freqs <= 0.15) hf_band = (freqs >= 0.15) & (freqs <= 0.4) # 4. Calculate Area Under Curve (AUC) lf_power = np.trapz(psd[lf_band], freqs[lf_band]) hf_power = np.trapz(psd[hf_band], freqs[hf_band]) return lf_power / hf_power # Usage # ratio = analyze_frequency_domain(my_intervals) # print(f"Burnout Index (LF/HF): {ratio:.2f}") COMMAND_BLOCK: from scipy.interpolate import interp1d from scipy.signal import welch def analyze_frequency_domain(nn_intervals): # 1. Resample: NN-intervals are irregularly spaced. We need 4Hz interpolation for FFT. time = np.cumsum(nn_intervals) / 1000.0 # seconds f_interp = interp1d(time, nn_intervals, kind='cubic') new_time = np.arange(time[0], time[-1], 0.25) # 4Hz nn_resampled = f_interp(new_time) # 2. Power Spectral Density using Welch's method freqs, psd = welch(nn_resampled, fs=4.0, nperseg=256) # 3. Define bands lf_band = (freqs >= 0.04) & (freqs <= 0.15) hf_band = (freqs >= 0.15) & (freqs <= 0.4) # 4. Calculate Area Under Curve (AUC) lf_power = np.trapz(psd[lf_band], freqs[lf_band]) hf_power = np.trapz(psd[hf_band], freqs[hf_band]) return lf_power / hf_power # Usage # ratio = analyze_frequency_domain(my_intervals) # print(f"Burnout Index (LF/HF): {ratio:.2f}") - LF (Low Frequency: 0.04 - 0.15 Hz): Reflects both Sympathetic (Fight or Flight) and Parasympathetic activity. - HF (High Frequency: 0.15 - 0.40 Hz): Reflects Parasympathetic (Rest and Digest) activity.