heart rate improvements

changelog.md

@@ -5,14 +5,18 @@
 - Added parameters that appear when attempting to generate results from the Participant and Intra-Group Functional Connectivity viewers and removed the non-functional placeholder parameters
 - Modified the Participant and Intra-Group Functional Connectivity analysis options to better perform their expected tasks. This remains as a BETA feature
 - Removed the existing Intra-Group Functional Connectivity option and replaced it with two new ones: Beta-Series Correlation and Spectral Coherence (epochs)
-- Updated the names of the methods provided for the Participant Functional Coneectivity Viewer to better match the actions they perform
+- Updated the names of the methods provided for the Participant Functional Connectivity Viewer to better match the actions they perform
 - Updated the warnings for the Functional Connectivity Viewers to better represent the challenges these analysis options now face
 - Added basic unit testing to hopefully prevent any accidental processing changes from occurring in the future
+- Added a new parameter to the PSP section on the right side of the screen: PSP_USE_HEART_RATE_BAND. This functions similarly to the existing SCI_USE_HEART_RATE_BAND
 - Added description text to the Inter-Group and Intra-Group Brain and Image Viewers, as well as the Functional Connectivity windows to explain what output can be expected
 - Removed image index 1 (Significance) from the Intra-Group Brain and Image Viewer as it is now provided more in depth with the Stats viewers
 - Modified the timeout when waiting for the application to close while performing updates down to a reasonable number
+- Modified the help messages for parameters in the SCI and PSP areas to better reflect how the parameters are used
 - Modified the heart rate calculation to be more precise and correct when dealing with good data and not messy data
 - Modified the heart rate calculation to not take only one channel in the data to use, but rather an average of channels. This still prefers short channels if they are present
+- Fixed an issue where performing PSP could ignore the bad channels that were marked by SCI and SNR
+- Fixed an issue where the heart rate calculation would disregard all calculations and fall back to an extremely rudimentary calculation
 - Fixed an issue that could prevent log file generation while the application was in the middle of an update
 - Fixed an issue where a rare crash could occur while the application was in the middle of an update
 - Fixed an issue that could have passed multiple conditions when generating an Intra-Group Stats image

flares.py

@@ -56,6 +56,11 @@ from scipy.spatial.distance import cdist
 from scipy.signal import welch, butter, filtfilt, periodogram # type: ignore
 from scipy.stats import pearsonr, zscore, ttest_1samp, ttest_ind, sem, t as t_dist # type: ignore
 
+from openpyxl import Workbook
+from openpyxl.styles import Font, PatternFill, Alignment
+from openpyxl.formatting.rule import ColorScaleRule
+from openpyxl.utils import get_column_letter
+
 import pywt # type: ignore
 import neurokit2 as nk # type: ignore
 
@@ -106,7 +111,7 @@ from src.shared.shareddata import PLATFORM_NAME, resource_path
 
 PRIMARY_COLORS = {
     "SCI only": "skyblue",      # Scalp Coupling Index (Standard MNE)
-    "SNR only": "lightgreen",   # Signal-to-Noise Ratio (Your original)
+    "SNR only": "lightgreen",   # Signal-to-Noise Ratio (Original)
     "PSP only": "salmon",       # Power Spectral Peak (Original Noise check)
     "Coeff_var only": "yellow", # Relative Noise (The coeff_var-only check)
     "Range only": "coral",      # Z-Swing (The Range Outlier check)
@@ -122,6 +127,42 @@ def get_category_color(label):
     """Returns the primary color if it's a single failure, otherwise gray."""
     return PRIMARY_COLORS.get(label, COMBINATION_COLOR)
 
+
+# direction: True = lower value is better (green), False = higher is better (green)
+QC_METRIC_DIRECTIONS = {
+    "n_bad_sci": True,
+    "n_bad_snr": True,
+    "n_bad_psp": True,
+    "n_bad_coeff_var": True,
+    "n_bad_mad": True,
+    "n_bad_psd_noise": True,
+    "n_bad_dropout": True,
+    "n_bad_channels_total": True,
+    "pct_bad_channels": True,
+    "total_processing_seconds": True,
+    "n_epochs_final": False,
+}
+
+# metrics with no inherent "good direction" - colored by deviation from the
+# group median instead (outliers flagged, not high/low values per se)
+QC_METRIC_DEVIATION_BASED = {"final_hr_bpm"}
+
+QC_METRIC_LABELS = {
+    "n_bad_sci": "Bad Channels - SCI",
+    "n_bad_snr": "Bad Channels - SNR",
+    "n_bad_psp": "Bad Channels - PSP",
+    "n_bad_coeff_var": "Bad Channels - Coeff. Var",
+    "n_bad_mad": "Bad Channels - MAD",
+    "n_bad_psd_noise": "Bad Channels - PSD Noise",
+    "n_bad_dropout": "Bad Channels - Dropout",
+    "n_bad_channels_total": "Bad Channels - Total (union)",
+    "pct_bad_channels": "% Channels Bad",
+    "final_hr_bpm": "Final Heart Rate (BPM)",
+    "n_epochs_final": "Epochs Retained",
+    "total_processing_seconds": "Processing Time (s)",
+}
+
+
 DOWNSAMPLE: bool
 DOWNSAMPLE_FREQUENCY: int
 
@@ -162,10 +203,11 @@ SNR_NOISE_LOW_FREQ: float
 SNR_NOISE_HIGH_FREQ: float
 
 PSP: bool
-PSP_TIME_WINDOW: int
-PSP_THRESHOLD: float
+PSP_USE_HEART_RATE_BAND: bool
 PSP_LOW_FREQ: float
 PSP_HIGH_FREQ: float
+PSP_TIME_WINDOW: int
+PSP_THRESHOLD: float
 
 COEFF_VAR: bool
 COEFF_VAR_THRESHOLD: int
@@ -260,6 +302,7 @@ GROUP: str = "Default"
 
 FOLDING_BYP: bool = False
 
+FEATURE_1: bool = False
 
 # Ensure that we are working in the directory of this file
 script_dir = os.path.dirname(os.path.abspath(__file__))
@@ -416,7 +459,7 @@ def process_participant_worker(file_path, file_params, file_metadata, result_que
 
 
 def process_multiple_participants(file_paths, file_params, file_metadata,
-                                   progress_queue=None, gui_queue=None, max_workers=6):
+                                   progress_queue=None, gui_queue=None, max_workers=6, qc_summary_path: str | None = None):
     ctx = mp.get_context("spawn")
     result_queue = ctx.Queue()
 
@@ -433,6 +476,8 @@ def process_multiple_participants(file_paths, file_params, file_metadata,
     duration_total = {"value": 0.0}
     success_count = {"value": 0}
     failed_stages = {"value": []}
+    qc_rows: list[dict[str, Any]] = []
+    qc_summary_path="qc_summary.xlsx"
 
     def elapsed_heartbeat():
         # Ticks once a second so the GUI can show a live-updating timer,
@@ -471,8 +516,20 @@ def process_multiple_participants(file_paths, file_params, file_metadata,
             if error is None:
                 duration_total["value"] += duration
                 success_count["value"] += 1
+                qc = result[-2] if isinstance(result, tuple) else None # accessing by index of return
+                if isinstance(qc, dict):
+                    qc["status"] = "success"
+                    qc["duration_seconds"] = round(duration, 2)
+                    qc_rows.append(qc)
             else:
                 failed_stages["value"].append(stage)
+                qc_rows.append({
+                    "file_path": res_path,
+                    "status": "FAILED",
+                    "error": error.splitlines()[0] if error else "unknown error",
+                    "duration_seconds": round(duration, 2),
+                    "total_processing_seconds": None,
+                })
             if gui_queue:
                 try:
                     gui_queue.put({
@@ -553,6 +610,14 @@ def process_multiple_participants(file_paths, file_params, file_metadata,
     process_multiple_participants._success_count = success_count["value"]
     process_multiple_participants._failed_stages = failed_stages["value"]
 
+    if qc_summary_path and qc_rows:
+        if FEATURE_1:
+            try:
+                write_qc_excel_summary(qc_rows, qc_summary_path)
+                logger.info(f"QC summary written to {qc_summary_path} ({len(qc_rows)} participant(s))")
+            except Exception as e:
+                logger.error(f"Failed to write QC summary: {e}")
+
     return results_by_file
 
 
@@ -1172,7 +1237,10 @@ def calculate_peak_power(data: BaseRaw, time_window: int = 3, threshold: float =
 
     # Identify channels that don't meet the provided threshold
     psp = scores.mean(axis=1)
-    data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), psp < threshold))
+    bad_channels = list(compress(cast(list[str], data.ch_names), psp < threshold))
+
+    existing_bads = set(data.info.get("bads", []))
+    data.info["bads"] = list(existing_bads | set(bad_channels))
 
     # Determine the colors based on the threshold, and create the figures
     color_stops = ([0.0, threshold, threshold+0.1, threshold+0.2, 1.0], [0.0, threshold, threshold, 1.0])
@@ -4409,17 +4477,22 @@ def _select_hr_source_channels(
     detection to already have succeeded.
     """
     n_channels = channel_data.shape[0]
-    scores = np.zeros(n_channels)
     nperseg = min(channel_data.shape[1], 2048)
     if nperseg < 8:
-        return scores  # too little data for a meaningful PSD
+        return np.zeros(n_channels)
 
-    for ch in range(n_channels):
-        freqs, psd = welch(channel_data[ch], fs=sfreq, nperseg=nperseg)
-        band_mask = (freqs >= cardiac_band[0]) & (freqs <= cardiac_band[1])
-        total_power = np.sum(psd)
-        scores[ch] = np.sum(psd[band_mask]) / total_power if total_power > 0 else 0.0
+    # single vectorized call across all channels, instead of one welch()
+    # call per channel - same math, same result, no Python-level loop
+    freqs, psd = welch(channel_data, fs=sfreq, nperseg=nperseg, axis=1)  # psd shape: (n_channels, n_freqs)
 
+    band_mask = (freqs >= cardiac_band[0]) & (freqs <= cardiac_band[1])
+    total_power = np.sum(psd, axis=1)
+    band_power = np.sum(psd[:, band_mask], axis=1)
+
+    scores = np.divide(
+        band_power, total_power,
+        out=np.zeros(n_channels), where=total_power > 0
+    )
     return scores
 
 
@@ -4525,6 +4598,104 @@ def short_channel_processing_for_hr(
     return sfreq, signal_trimmed, times_trimmed
 
 
+def reconcile_heart_rate_estimates(
+    mean_hr_scipy: float,
+    psd_confidence: float,
+    mean_hr_nk: float,
+    mode_hr_nk: float,
+    agreement_tolerance_bpm: float = 10.0,
+    psd_confidence_threshold: float = 3.0,
+) -> tuple[float, bool, str]:
+    """
+    Combines three independent heart rate estimates - PSD spectral peak
+    (scipy), NeuroKit's cleaned/interpolated mean, and NeuroKit's mode -
+    into one final value, using majority agreement rather than a fixed
+    pairwise override rule.
+
+    Logic, in order:
+    1. If all three estimates agree within agreement_tolerance_bpm of each
+       other, average them - strongest possible evidence, no single method
+       is being trusted over the others.
+    2. Otherwise, check if any TWO of the three agree with each other -
+       if so, average that agreeing pair and discard the outlier. Two
+       independent methods landing on the same value by coincidence is
+       unlikely; the third is more likely the one that's wrong.
+    3. If no two agree at all, fall back to whichever single estimate is
+       most trustworthy: the PSD estimate if its confidence clears
+       psd_confidence_threshold (a genuinely sharp, unambiguous spectral
+       peak), otherwise the NeuroKit mode (more robust than its mean to
+       a residual minority of bad samples, per _mode_hr's reasoning).
+
+    Returns
+    -------
+    tuple[float, bool, str]
+        - float: final reconciled heart rate (BPM).
+        - bool: True if any disagreement/overruling occurred (for plotting/logging).
+        - str: human-readable explanation of which path was taken, for logs.
+    """
+    estimates = {
+        "psd": mean_hr_scipy,
+        "nk_mean": mean_hr_nk,
+        "nk_mode": mode_hr_nk,
+    }
+
+    def _fmt(d: dict[str, float]) -> str:
+        return ", ".join(f"{k}={v:.1f}" for k, v in d.items())
+    
+    pairs = [("psd", "nk_mean"), ("psd", "nk_mode"), ("nk_mean", "nk_mode")]
+    agreeing_pairs = [
+        (a, b) for a, b in pairs
+        if abs(estimates[a] - estimates[b]) <= agreement_tolerance_bpm
+    ]
+
+    if len(agreeing_pairs) == 3:
+        final = float(np.mean(list(estimates.values())))
+        return final, False, f"All three estimates agree (within {agreement_tolerance_bpm} BPM) - averaged: {_fmt(estimates)}"
+
+    if len(agreeing_pairs) >= 1:
+        a, b = agreeing_pairs[0]
+        final = float((estimates[a] + estimates[b]) / 2.0)
+        outlier = [k for k in estimates if k not in (a, b)][0]
+        return final, True, (
+            f"{a} and {b} agree ({estimates[a]:.1f}, {estimates[b]:.1f}); "
+            f"{outlier} is an outlier ({estimates[outlier]:.1f}) - discarded."
+        )
+
+    # No two estimates agree at all - fall back to the single most trustworthy one
+    if psd_confidence >= psd_confidence_threshold:
+        return mean_hr_scipy, True, (
+            f"No two estimates agree ({_fmt(estimates)}); PSD peak is clear "
+            f"(confidence={psd_confidence:.2f}) - trusting PSD alone."
+        )
+    else:
+        return mode_hr_nk, True, (
+            f"No two estimates agree ({_fmt(estimates)}); PSD peak is ambiguous "
+            f"(confidence={psd_confidence:.2f} < {psd_confidence_threshold}) - "
+            f"trusting NeuroKit mode alone (more robust to residual dips than its mean)."
+        )
+
+    
+def _mode_hr(hr_clean: NDArray[float64], bin_width_bpm: float = 2.0) -> float:
+    """
+    Histogram-based mode of the HR trace: the center of the most frequently
+    occurring bin. More robust than a plain mean to a minority of corrupted
+    (missed-beat) dips, since those dips only need to avoid being the
+    single largest cluster - unlike a mean, which every dip pulls down
+    proportionally regardless of how rare it is.
+
+    bin_width_bpm : float, default 2.0
+        Histogram bin width. Too narrow and there's no meaningful mode
+        (every value nearly unique); too wide and you lose real precision
+        in the estimate. 2 BPM is a reasonable starting point - worth
+        checking against your actual HR distributions.
+    """
+    if len(hr_clean) == 0:
+        return float('nan')
+    bins = np.arange(hr_clean.min(), hr_clean.max() + bin_width_bpm, bin_width_bpm)
+    counts, edges = np.histogram(hr_clean, bins=bins)
+    mode_bin_idx = np.argmax(counts)
+    return float((edges[mode_bin_idx] + edges[mode_bin_idx + 1]) / 2.0)
+
 
 def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64], hr_low_freq, hr_high_freq, max_low_hr, max_high_hr, smoothing_window_hr, short_channels) -> tuple[NDArray[float64], float]:
     """
@@ -4575,12 +4746,14 @@ def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64]
     smoothed_values[spikes] = smoothed_spikes[spikes]
     hr_smooth_nk = cast(NDArray[float64], smoothed_values.to_numpy())
     mean_hr_nk = hr_smooth_nk.mean()
-
+    mode_hr_nk = _mode_hr(hr_clean, bin_width_bpm=2.0)
     logger.info("Original HR min/max: %f, %f", hr_clean.min(), hr_clean.max())
     logger.info("Smoothed HR min/max:%f, %f", hr_smooth_nk.min(), hr_smooth_nk.max())
-    logger.info(f"Estimated mean HR nk: {mean_hr_nk:.1f} BPM")
+    logger.info(f"Estimated mean HR nk: {mean_hr_nk:.1f} BPM, mode HR nk: {mode_hr_nk:.1f} BPM")
+
+
+    return hr_smooth_nk, mean_hr_nk, mode_hr_nk
 
-    return hr_smooth_nk, mean_hr_nk
 
 
 def calculate_heart_rate_scipy(
@@ -4678,8 +4851,11 @@ def plot_heart_rate(
     mean_hr_scipy: float,
     hr_smooth_nk: NDArray[floating[Any]],
     mean_hr_nk: float,
+    mode_hr_nk: float,
+    final_hr: float,
     times_trimmed: NDArray[floating[Any]],
     overruled: bool,
+    reconciliation_note: str,
     hr_window: int
 ) -> tuple[Figure, Figure]:
     """
@@ -4694,15 +4870,26 @@ def plot_heart_rate(
     freq_range_scipy : np.ndarray[Any, np.dtype[np.bool_]]
         Boolean mask indicating the heart rate frequency range used in PSD.
     mean_hr_scipy : float
-        Mean heart rate estimated from SciPy PSD peak.
+        Heart rate estimated from SciPy PSD peak (cluster centroid).
     hr_smooth_nk : NDArray[floating[Any]]
         Smoothed instantaneous heart rate from NeuroKit2.
     mean_hr_nk : float
-        Mean heart rate estimated from NeuroKit2 data.
+        Mean heart rate estimated from NeuroKit2's cleaned time series.
+    mode_hr_nk : float
+        Mode (most frequent value) of NeuroKit2's heart rate distribution -
+        more robust than the mean to a residual minority of missed-beat dips.
+    final_hr : float
+        The final reconciled heart rate value, combining all three estimates
+        (see reconcile_heart_rate_estimates).
     times_trimmed : NDArray[floating[Any]]
         Time points corresponding to hr_smooth_nk values.
-    overruled: bool
-        True if the heart rate from NeuroKit2 is overriding the results from the PSD.
+    overruled : bool
+        True if the final reconciled value differs from a simple average of
+        all three (i.e. an outlier was discarded or a single estimate was
+        trusted alone).
+    reconciliation_note : str
+        Human-readable explanation of which reconciliation path was taken -
+        shown directly on the plot rather than a generic "was overruled" note.
 
     Returns
     -------
@@ -4724,32 +4911,26 @@ def plot_heart_rate(
     ax1.grid(True) # type: ignore
 
     # Was the value we reported here correct for the data on the graph or was it overruled?
-    if overruled:
-        note = (
-            '\n'
-            'Note: Calculation was bad!\n'
-            'Data has been set to match\n'
-            'the value from NeuroKit2.'
-        )
-        phantom = Line2D([0], [0], color='none', label=note)
-        handles, _ = ax1.get_legend_handles_labels()
-        ax1.legend(handles=handles + [phantom]) # type: ignore
-
-    else:
-        ax1.legend() # type: ignore
+    note_prefix = 'Reconciliation (outlier discarded):' if overruled else 'Reconciliation (estimates agreed):'
+    note = f"\n{note_prefix}\n{reconciliation_note}"
+    phantom = Line2D([0], [0], color='none', label=note)
+    handles, _ = ax1.get_legend_handles_labels()
+    ax1.legend(handles=handles + [phantom], fontsize=8)
     plt.close(fig1)
 
     # Create the second plot showing the rolling heart rate, as well as the two averages that were calculated
     logger.info("Creating the figure...")
     fig2, ax2 = plt.subplots(figsize=(14, 6)) # type: ignore
-    ax2.plot(times_trimmed, hr_smooth_nk, label='Instantaneous HR (NeuroKit2)', color='blue', alpha=0.7) # type: ignore
-    ax2.axhline(mean_hr_nk, color='red', linestyle='--', label=f'Mean HR NeuroKit2: {mean_hr_nk:.1f} BPM') # type: ignore
-    ax2.axhline(mean_hr_scipy, color='orange', linestyle=':', label=f'SciPy Welch PSD (HP filtered): {mean_hr_scipy:.1f} BPM') # type: ignore
-    ax2.set_xlabel('Time (seconds)') # type: ignore
-    ax2.set_ylabel('Heart Rate (BPM)') # type: ignore
-    ax2.set_title('Heart Rate Estimates Comparison') # type: ignore
-    ax2.legend() # type: ignore
-    ax2.grid(True) # type: ignore
+    ax2.plot(times_trimmed, hr_smooth_nk, label='Instantaneous HR (NeuroKit2)', color='blue', alpha=0.7)
+    ax2.axhline(mean_hr_nk, color='steelblue', linestyle='--', alpha=0.7, label=f'NK Mean: {mean_hr_nk:.1f} BPM')
+    ax2.axhline(mode_hr_nk, color='purple', linestyle='--', alpha=0.7, label=f'NK Mode: {mode_hr_nk:.1f} BPM')
+    ax2.axhline(mean_hr_scipy, color='orange', linestyle=':', label=f'PSD Estimate: {mean_hr_scipy:.1f} BPM')
+    ax2.axhline(final_hr, color='green', linestyle='-', linewidth=2, label=f'Final (Reconciled) HR: {final_hr:.1f} BPM')
+    ax2.set_xlabel('Time (seconds)')
+    ax2.set_ylabel('Heart Rate (BPM)')
+    ax2.set_title('Heart Rate Estimates Comparison')
+    ax2.legend(fontsize=9)
+    ax2.grid(True)
     fig2.tight_layout()
     plt.close(fig2)
 
@@ -4944,75 +5125,43 @@ def find_bad_channels_coeff_var(raw, coeff_var_threshold=25.0):
 
 
 
-def hr_calc(raw, seconds_to_strip_hr, l_freq, h_freq, search_min, search_max, max_low_hr, max_high_hr, smoothing_window_hr, hr_window, short_channels, short_channels_threshold, verbosity, psd_confidence_threshold: float = 3.0):
+def hr_calc(raw, seconds_to_strip_hr, l_freq, h_freq, search_min, search_max, max_low_hr, max_high_hr, smoothing_window_hr, hr_window, short_channels, short_channels_threshold, verbosity, psd_confidence_threshold: float = 3.0, band_halfwidth_hz: float = 0.3):
     if short_channels:
         short_chans = get_short_channels(raw, max_dist=short_channels_threshold)
     else:
         short_chans = None
     sfreq, signal_trimmed, times_trimmed = short_channel_processing_for_hr(raw, short_chans, seconds_to_strip_hr=seconds_to_strip_hr, verbosity=verbosity)
-    hr_smooth_nk, mean_hr_nk = calculate_heart_rate_neurokit(sfreq, signal_trimmed, hr_low_freq=l_freq, hr_high_freq=h_freq, max_low_hr=max_low_hr, max_high_hr=max_high_hr, smoothing_window_hr=smoothing_window_hr, short_channels=short_channels)
+    hr_smooth_nk, mean_hr_nk, mode_hr_nk = calculate_heart_rate_neurokit(sfreq, signal_trimmed, hr_low_freq=l_freq, hr_high_freq=h_freq, max_low_hr=max_low_hr, max_high_hr=max_high_hr, smoothing_window_hr=smoothing_window_hr, short_channels=short_channels)
     freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy, psd_confidence = calculate_heart_rate_scipy(sfreq, signal_trimmed, search_min=search_min, search_max=search_max)
-    
-    overruled = False
-    disagreement = abs(mean_hr_scipy - mean_hr_nk) > 15
 
-    if disagreement:
-        if psd_confidence >= psd_confidence_threshold:
-            logger.info(f"HR estimates disagree ({mean_hr_scipy:.1f} vs {mean_hr_nk:.1f} BPM) - "
-                        f"PSD peak is clear (confidence={psd_confidence:.2f}), trusting PSD. Overruling NeuroKit.")
-            mean_hr_nk = mean_hr_scipy
-            overruled = True
-        else:
-            logger.info(f"HR estimates disagree ({mean_hr_scipy:.1f} vs {mean_hr_nk:.1f} BPM) - "
-                        f"PSD peak is ambiguous (confidence={psd_confidence:.2f} < {psd_confidence_threshold}), "
-                        f"trusting NeuroKit instead. Overruling PSD.")
-            mean_hr_scipy = mean_hr_nk
-            overruled = True
+    final_hr, overruled, reconciliation_note = reconcile_heart_rate_estimates(
+        mean_hr_scipy, psd_confidence, mean_hr_nk, mode_hr_nk, agreement_tolerance_bpm=hr_window
+    )
+    logger.info(f"HR reconciliation: {reconciliation_note}")
+    logger.info(f"Final heart rate: {final_hr:.1f} BPM")
 
-    hr1, hr2 = plot_heart_rate(freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy, hr_smooth_nk, mean_hr_nk, times_trimmed, overruled, hr_window=hr_window)
+    mean_hr_scipy = final_hr
+
+    hr1, hr2 = plot_heart_rate(
+        freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy,
+        hr_smooth_nk, mean_hr_nk, mode_hr_nk, final_hr,
+        times_trimmed, overruled, reconciliation_note, hr_window=hr_window
+    )
    
     fig = raw.compute_psd().plot(show=False)
-    raw_filtered = raw.copy().filter(0.5, 3, fir_design='firwin')
-    sfreq = raw.info['sfreq']
-    data = raw_filtered.get_data()
-    channel_names = raw.ch_names
 
-    # --- Parameters for PSD ---
-    desired_bin_hz = 0.1
-    nperseg = int(sfreq / desired_bin_hz)
-    hr_range = (search_min, search_max)
+    # Targeted frequency band for downstream calculations (SCI, PSP, etc.),
+    # derived directly from the reconciled heart rate estimate above -
+    # replaces a previous, disconnected all-channel median calculation that
+    # never used the short-channel selection or reconciliation logic at all.
+    hr_freq = final_hr / 60.0  # BPM -> Hz
+    low = hr_freq - band_halfwidth_hz
+    high = hr_freq + band_halfwidth_hz
 
-    f, Pxx = welch(data, fs=sfreq, nperseg=nperseg, axis=1)  # data: (n_channels, n_samples)
-    mask = (f >= hr_range[0] / 60) & (f <= hr_range[1] / 60)
-    f_masked = f[mask]
-    Pxx_masked = Pxx[:, mask]  # (n_channels, n_freq_in_range)
- 
-    hr_all_channels = np.full(Pxx_masked.shape[0], np.nan)
-    if Pxx_masked.shape[1] >= 3:
-        # same "strictly greater than both neighbors" local-max definition as
-        # the original per-channel loop, vectorized across all channels at once
-        interior = Pxx_masked[:, 1:-1]
-        left = Pxx_masked[:, :-2]
-        right = Pxx_masked[:, 2:]
-        is_local_peak = (interior > left) & (interior > right)
- 
-        for ch in range(Pxx_masked.shape[0]):
-            peak_offsets = np.where(is_local_peak[ch])[0]
-            if len(peak_offsets) == 0:
-                continue
-            candidate_idx = peak_offsets + 1  # shift back into Pxx_masked indexing
-            best_idx = candidate_idx[np.argmax(Pxx_masked[ch, candidate_idx])]
-            hr_all_channels[ch] = f_masked[best_idx] * 60  # bpm
- 
-    hr_all_channels = hr_all_channels[~np.isnan(hr_all_channels)]
-    hr_mode = np.round(np.median(hr_all_channels))  # Use median if some NaNs
+    logger.info(f"SCI/PSP target band: {low:.3f}-{high:.3f} Hz "
+                f"({final_hr:.1f} +/- {band_halfwidth_hz*60:.0f} BPM)")
 
-    print(f"Estimated Heart Rate: {hr_mode} bpm")
-
-    hr_freq = hr_mode / 60  # Hz
-    low = hr_freq - 0.3
-    high = hr_freq + 0.3
-    return fig, hr1, hr2, low, high
+    return fig, hr1, hr2, low, high, final_hr
 
 
 
@@ -5353,7 +5502,7 @@ def _png_worker(png_queue: Queue, fig_bytes_dict: dict, dpi: int = 100):
 
 
 
-def initial_setup():
+def initial_setup(file_path):
     timings = {}
     step_start = time.perf_counter()
     config_dict = {
@@ -5363,20 +5512,21 @@ def initial_setup():
     }
 
     fig_bytes_dict: dict[str, bytes] = {}
+    qc: dict[str, Any] = {"file_path": file_path}
     png_queue: Queue = Queue()
     png_thread = threading.Thread(
         target=_png_worker, args=(png_queue, fig_bytes_dict), daemon=True
     )
     png_thread.start()
 
-    return fig_bytes_dict, config_dict, png_queue, timings, step_start
+    return fig_bytes_dict, config_dict, png_queue, timings, step_start, qc
 
 
 
 def _enqueue(label, fig, png_queue):
     if fig is None:
         return
-    plt.close(fig)  # detach from pyplot's global registry - main thread only
+    plt.close(fig)
     png_queue.put((label, fig))
 
 
@@ -5392,11 +5542,12 @@ def process_participant(file_path, file_start, progress_callback=None):
 
     print(f"File was started with {time.time() - file_start:2f} seconds elapsed.")
     # Step 0: Setting up
-    fig_bytes_dict, config_dict, png_queue, timings, step_start = initial_setup()
+    fig_bytes_dict, config_dict, png_queue, timings, step_start, qc = initial_setup(file_path)
     step_start = lap(step_start, timings, "Step 0")
 
     # Step 1: Preprocessing
     raw = load_snirf(file_path=file_path, downsample_frequency=DOWNSAMPLE_FREQUENCY, verbosity=VERBOSITY)
+    qc["n_channels_loaded"] = raw.info['nchan']
     fig_raw = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Loaded Raw", show=False)
     _enqueue("Loaded Raw Data", fig_raw, png_queue)
     if progress_callback: progress_callback(1)
@@ -5432,8 +5583,9 @@ def process_participant(file_path, file_start, progress_callback=None):
     step_start = lap(step_start, timings, "Step 4")
 
     # Step 5: Heart Rate
+    qc["heart_rate_ran"] = HEART_RATE and not FOLDING_BYP
     if HEART_RATE and not FOLDING_BYP:
-        fig, hr1, hr2, low, high = hr_calc(
+        fig, hr1, hr2, low, high, final_hr = hr_calc(
             raw, 
             seconds_to_strip_hr=SECONDS_TO_STRIP_HR, 
             l_freq=HR_LOW_FREQ, 
@@ -5448,6 +5600,7 @@ def process_participant(file_path, file_start, progress_callback=None):
             short_channels_threshold=SHORT_CHANNELS_THRESHOLD,
             verbosity=VERBOSITY
         )
+        qc["final_hr_bpm"] = round(final_hr, 1)
         _enqueue("Power Spectral Density", fig, png_queue)
         _enqueue('Heart Rate - PSD', hr1, png_queue)
         _enqueue('Heart Rate - Time', hr2, png_queue)
@@ -5464,6 +5617,7 @@ def process_participant(file_path, file_start, progress_callback=None):
             bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, l_freq=SCI_LOW_FREQ, h_freq=SCI_HIGH_FREQ, time_window=SCI_TIME_WINDOW, threshold=SCI_THRESHOLD)
         _enqueue("Scalp Coupling Index Heatmap", fig_sci_1, png_queue)
         _enqueue("Scalp Coupling Index Binary Heatmap", fig_sci_2, png_queue)
+    qc["n_bad_sci"] = len(bad_sci)
     if progress_callback: progress_callback(6)
     logger.info("Step 6 Completed.")
     step_start = lap(step_start, timings, "Step 6")
@@ -5473,6 +5627,7 @@ def process_participant(file_path, file_start, progress_callback=None):
     if SNR and not FOLDING_BYP:
         bad_snr, fig_snr = calculate_signal_noise_ratio(raw)
         _enqueue("Signal To Noise Ratio", fig_snr, png_queue)
+    qc["n_bad_snr"] = len(bad_snr)
     if progress_callback: progress_callback(7)
     logger.info("Step 7 Completed.")
     step_start = lap(step_start, timings, "Step 7")
@@ -5480,9 +5635,13 @@ def process_participant(file_path, file_start, progress_callback=None):
     # Step 8: Peak Spectral Power
     bad_psp = []
     if PSP and not FOLDING_BYP:
-        bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=PSP_LOW_FREQ, h_freq=PSP_HIGH_FREQ)
+        if HEART_RATE and PSP_USE_HEART_RATE_BAND:
+            bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=low, h_freq=high)
+        else:
+            bad_psp, fig_psp1, fig_psp2 = calculate_peak_power(raw, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=PSP_LOW_FREQ, h_freq=PSP_HIGH_FREQ)
         _enqueue("Peak Spectral Power Heatmap", fig_psp1, png_queue)
         _enqueue("Peak Spectral Power Binary Heatmap", fig_psp2, png_queue)
+    qc["n_bad_psp"] = len(bad_psp)
     if progress_callback: progress_callback(8)
     logger.info("Step 8 Completed.")
     step_start = lap(step_start, timings, "Step 8")
@@ -5492,6 +5651,7 @@ def process_participant(file_path, file_start, progress_callback=None):
     if COEFF_VAR and not FOLDING_BYP:
         bad_coeff_var, fig_coeff_var = find_bad_channels_coeff_var(raw, coeff_var_threshold=COEFF_VAR_THRESHOLD)
         _enqueue('Coefficient of Variation', fig_coeff_var, png_queue)
+    qc["n_bad_coeff_var"] = len(bad_coeff_var)
     if progress_callback: progress_callback(9)
     logger.info("Step 9 Completed.")
     step_start = lap(step_start, timings, "Step 9")
@@ -5501,6 +5661,7 @@ def process_participant(file_path, file_start, progress_callback=None):
     if MAD and not FOLDING_BYP:
         bad_amplitude_range, fig_range = find_bad_channels_by_amplitude_range(raw, threshold=MAD_THRESHOLD)
         _enqueue('Median Absolute Deviation', fig_range, png_queue)
+    qc["n_bad_mad"] = len(bad_amplitude_range)
     if progress_callback: progress_callback(10)
     logger.info("Step 10 Completed.")
     step_start = lap(step_start, timings, "Step 10")
@@ -5510,6 +5671,7 @@ def process_participant(file_path, file_start, progress_callback=None):
     if PSD_NOISE and not FOLDING_BYP:
         bad_noise, fig_noise = detect_spectral_noise_spike(raw, db_limit=DB_LIMIT, freq_div=TARGET_FREQ_DIV, min_freq=PSD_MIN_FREQ, target_bandwith=PSD_TARGET_BANDWIDTH)
         _enqueue('Power Spectral Density Noise', fig_noise, png_queue)
+    qc["n_bad_psd_noise"] = len(bad_noise)
     if progress_callback: progress_callback(11)
     logger.info("Step 11 Completed.")
     step_start = lap(step_start, timings, "Step 11")
@@ -5519,13 +5681,20 @@ def process_participant(file_path, file_start, progress_callback=None):
     if SENSOR_DROPOUT and not FOLDING_BYP:
         bad_disp, fig_disp = detect_sensor_dropout(raw, threshold_ratio=SENSOR_DROPOUT_VARIANCE_THRESHOLD)
         _enqueue('Sensor Dropout', fig_disp, png_queue)
+    qc["n_bad_dropout"] = len(bad_disp)
     if progress_callback: progress_callback(12)
     logger.info("Step 12 Completed.")
     step_start = lap(step_start, timings, "Step 12")
 
     # Step 13: Bad Channels Handling
+    qc["n_bad_channels_total"] = 0
+    qc["pct_bad_channels"] = 0.0
+    qc["bad_channels_handling"] = "None"
     if BAD_CHANNELS_HANDLING != "None" and not FOLDING_BYP:
         raw, fig_dropped, fig_raw_before, bad_channels = mark_bads(raw, bad_sci, bad_snr, bad_psp, bad_coeff_var, bad_amplitude_range, bad_noise, bad_disp)
+        qc["n_bad_channels_total"] = len(bad_channels)
+        qc["pct_bad_channels"] = round(100 * len(bad_channels) / qc["n_channels_loaded"], 1) if qc["n_channels_loaded"] else 0.0
+        qc["bad_channels_handling"] = BAD_CHANNELS_HANDLING
         if fig_dropped and fig_raw_before is not None:
             _enqueue("Bad Channels by Method", fig_dropped, png_queue)
             _enqueue("Bad Channels Data", fig_raw_before, png_queue)
@@ -5625,6 +5794,8 @@ def process_participant(file_path, file_start, progress_callback=None):
     step_start = lap(step_start, timings, "Step 20")
 
     # Step 21: Epoch Calculations
+    epochs = None
+    qc["n_epochs_final"] = 0
     if EPOCHS and EVENTS and not FOLDING_BYP:
         epochs = epochs_calculations(
             raw_haemo_evnt, 
@@ -5639,6 +5810,7 @@ def process_participant(file_path, file_start, progress_callback=None):
             reject_hbo_threshold=dict(hbo=REJECT_HBO_THRESHOLD),
             png_queue=png_queue
         )
+        qc["n_epochs_final"] = len(epochs) if epochs is not None else 0
     if progress_callback: progress_callback(21)
     logger.info("Step 21 Completed.")
     step_start = lap(step_start, timings, "Step 21")
@@ -5714,12 +5886,13 @@ def process_participant(file_path, file_start, progress_callback=None):
     step_start = lap(step_start, timings, "Step 28")
 
     # Step 28.5: Return the results
+    qc["total_processing_seconds"] = round(sum(timings.values()), 2)
     logger.info("Step timings:")
     for name, elapsed in timings.items():
         logger.info(f"  {name:<25} {elapsed:7.3f}s")
 
     logger.info(f"Total processing time: {sum(timings.values()):.3f}s")
-    return raw_haemo, epochs, df_cha, df_roi, df_design_matrix, config_dict, fig_bytes_dict, contrast_results_dict, roi_channel_map, True
+    return raw_haemo, epochs, df_cha, df_roi, df_design_matrix, config_dict, fig_bytes_dict, contrast_results_dict, roi_channel_map, qc, True
 
 
 
@@ -6593,6 +6766,134 @@ def peak_power_fast(
     return raw, scores, times
 
 
+def write_qc_excel_summary(qc_rows: list[dict], output_path: str) -> None:
+    """
+    Writes one Excel workbook summarizing QC metrics across all participants
+    in a batch run - participants as columns, metrics as rows, each row
+    color-scaled green (best) to red (worst) with direction-aware coloring,
+    plus live summary formulas (mean/median/min/max/worst participant) per row.
+
+    Failed participants are shown in a separate, clearly labeled block so
+    they don't distort the color scale of successful participants' numbers.
+    """
+    successes = [r for r in qc_rows if r.get("status") == "success"]
+    failures = [r for r in qc_rows if r.get("status") == "FAILED"]
+
+    wb = Workbook()
+    ws = wb.active
+    ws.title = "QC Summary"
+
+    header_font = Font(name="Arial", bold=True, size=11)
+    label_font = Font(name="Arial", size=10)
+    body_font = Font(name="Arial", size=10)
+    fail_fill = PatternFill(start_color="FFC7CE", end_color="FFC7CE", fill_type="solid")
+
+    metric_keys = [k for k in QC_METRIC_LABELS if k in QC_METRIC_DIRECTIONS or k in QC_METRIC_DEVIATION_BASED]
+    n_participants = len(successes)
+
+    # --- Header row ---
+    ws.cell(row=1, column=1, value="Metric").font = header_font
+    for i, row in enumerate(successes):
+        col = i + 2
+        c = ws.cell(row=1, column=col, value=row.get("file_path", f"P{i+1}"))
+        c.font = header_font
+        c.alignment = Alignment(horizontal="center", wrap_text=True)
+
+    summary_start_col = n_participants + 3  # one blank column, then summary block
+    for j, label in enumerate(["Mean", "Median", "Min", "Max", "Worst Participant"]):
+        c = ws.cell(row=1, column=summary_start_col + j, value=label)
+        c.font = header_font
+        c.alignment = Alignment(horizontal="center", wrap_text=True)
+
+    # --- Metric rows ---
+    for r_idx, key in enumerate(metric_keys):
+        row_num = r_idx + 2
+        ws.cell(row=row_num, column=1, value=QC_METRIC_LABELS[key]).font = label_font
+
+        for i, row in enumerate(successes):
+            col = i + 2
+            val = row.get(key)
+            c = ws.cell(row=row_num, column=col, value=val)
+            c.font = body_font
+
+        if n_participants == 0:
+            continue
+
+        first_col_letter = get_column_letter(2)
+        last_col_letter = get_column_letter(n_participants + 1)
+        data_range = f"{first_col_letter}{row_num}:{last_col_letter}{row_num}"
+
+        mean_col, median_col, min_col, max_col, worst_col = (
+            summary_start_col, summary_start_col + 1, summary_start_col + 2,
+            summary_start_col + 3, summary_start_col + 4,
+        )
+
+        ws.cell(row=row_num, column=mean_col, value=f"=AVERAGE({data_range})").font = body_font
+        ws.cell(row=row_num, column=median_col, value=f"=MEDIAN({data_range})").font = body_font
+        ws.cell(row=row_num, column=min_col, value=f"=MIN({data_range})").font = body_font
+        ws.cell(row=row_num, column=max_col, value=f"=MAX({data_range})").font = body_font
+
+        header_range = f"{first_col_letter}1:{last_col_letter}1"
+        if key in QC_METRIC_DEVIATION_BASED:
+            # "worst" = furthest from the row's own median
+            worst_formula = (
+                f"=INDEX({header_range},MATCH(MAX(ABS({data_range}-MEDIAN({data_range}))),"
+                f"ABS({data_range}-MEDIAN({data_range})),0))"
+            )
+        elif QC_METRIC_DIRECTIONS[key]:  # lower is better -> worst = max
+            worst_formula = f"=INDEX({header_range},MATCH(MAX({data_range}),{data_range},0))"
+        else:  # higher is better -> worst = min
+            worst_formula = f"=INDEX({header_range},MATCH(MIN({data_range}),{data_range},0))"
+
+        ws.cell(row=row_num, column=worst_col, value=worst_formula).font = body_font
+
+        # --- Color scale, direction-aware ---
+        if key in QC_METRIC_DEVIATION_BASED:
+            # color by |value - row median| via a helper column pattern isn't
+            # natively supported by ColorScaleRule (it colors raw cell values,
+            # not a derived formula) - approximate by centering the 3-color
+            # scale on the row's own values, which still highlights the
+            # extremes/outliers visually even without true deviation coloring.
+            rule = ColorScaleRule(
+                start_type="min", start_color="63BE7B",
+                mid_type="percentile", mid_value=50, mid_color="FFEB84",
+                end_type="max", end_color="F8696B",
+            )
+        elif QC_METRIC_DIRECTIONS[key]:  # lower is better: green=min, red=max
+            rule = ColorScaleRule(
+                start_type="min", start_color="63BE7B",
+                mid_type="percentile", mid_value=50, mid_color="FFEB84",
+                end_type="max", end_color="F8696B",
+            )
+        else:  # higher is better: green=max, red=min
+            rule = ColorScaleRule(
+                start_type="min", start_color="F8696B",
+                mid_type="percentile", mid_value=50, mid_color="FFEB84",
+                end_type="max", end_color="63BE7B",
+            )
+        ws.conditional_formatting.add(data_range, rule)
+
+    # --- Failed participants block, separate and clearly marked ---
+    if failures:
+        fail_row_start = len(metric_keys) + 4
+        ws.cell(row=fail_row_start, column=1, value="FAILED PARTICIPANTS").font = Font(name="Arial", bold=True, size=12, color="CC0000")
+        for i, row in enumerate(failures):
+            r = fail_row_start + 1 + i
+            path_cell = ws.cell(row=r, column=1, value=row.get("file_path", "unknown"))
+            path_cell.font = body_font
+            path_cell.fill = fail_fill
+            err_cell = ws.cell(row=r, column=2, value=row.get("error", "unknown error"))
+            err_cell.font = body_font
+            err_cell.fill = fail_fill
+
+    ws.column_dimensions['A'].width = 32
+    for i in range(n_participants):
+        ws.column_dimensions[get_column_letter(i + 2)].width = 14
+    ws.freeze_panes = "B2"
+
+    wb.save(output_path)
+
+
 if __name__ == "__main__":
     print("This file has no functionality when not used in tandem with the FLARES application.")
         

main.py

@@ -129,10 +129,10 @@ SECTIONS = [
         "params": [
             {"name": "SCI", "default": True, "type": bool, "advanced": False, "help": "Calculate and mark channels bad based on their Scalp Coupling Index. This metric calculates the quality of the connection between the optode and the scalp."},
             {"name": "SCI_USE_HEART_RATE_BAND", "default": True, "type": bool, "depends_on": [{"parent_name": "SCI"}, {"parent_name": "HEART_RATE"}], "advanced": False, "help": "Adjust the SCI frequency band using the participant's estimated heart rate."},
-            {"name": "SCI_LOW_FREQ", "default": 0.7, "type": float, "depends_on": "SCI_USE_HEART_RATE_BAND", "depends_value": False,"advanced": True, "help": "Lower frequency cutoff for SCI bandpass filtering (Hz)."},
-            {"name": "SCI_HIGH_FREQ", "default": 1.5, "type": float, "depends_on": "SCI_USE_HEART_RATE_BAND", "depends_value": False, "advanced": True, "help": "Upper frequency cutoff for SCI bandpass filtering (Hz)."},
+            {"name": "SCI_LOW_FREQ", "default": 0.7, "type": float, "depends_on": "SCI_USE_HEART_RATE_BAND", "depends_value": False,"advanced": True, "help": "Lower frequency bound for the signal band used in SCI calculation (Hz)."},
+            {"name": "SCI_HIGH_FREQ", "default": 1.5, "type": float, "depends_on": "SCI_USE_HEART_RATE_BAND", "depends_value": False, "advanced": True, "help": "Lower frequency bound for the signal band used in SCI calculation (Hz)."},
             {"name": "SCI_TIME_WINDOW", "default": 3, "type": int, "depends_on": "SCI", "advanced": False, "help": "Duration of each independent SCI calculation window in seconds."},
-            {"name": "SCI_THRESHOLD", "default": 0.6, "type": float, "depends_on": "SCI", "advanced": False, "help": "SCI threshold on a scale of 0-1. Channels below this value are marked bad."},
+            {"name": "SCI_THRESHOLD", "default": 0.6, "type": float, "depends_on": "SCI", "advanced": False, "help": "SCI threshold on a scale of 0-1. Channels below this value will be marked bad."},
         ]
     },
     {
@@ -150,10 +150,11 @@ SECTIONS = [
         "title": "Peak Spectral Power",
         "params": [
             {"name": "PSP", "default": True, "type": bool, "advanced": False, "help": "Calculate and mark channels bad based on their Peak Spectral Power. This metric calculates the amplitude or strength of the most prominent frequency component in a specified spectral range."},
-            {"name": "PSP_TIME_WINDOW", "default": 3, "type": int, "depends_on": "PSP", "advanced": False, "help": "Length of each independent PSP calculation window in seconds."},
-            {"name": "PSP_THRESHOLD", "default": 0.1, "type": float, "depends_on": "PSP", "advanced": False, "help": "Channels with average PSP values below this threshold will be marked as bad."},
-            {"name": "PSP_LOW_FREQ", "default": 0.7, "type": float, "depends_on": "PSP", "advanced": True, "help": "Lower frequency cutoff for PSP bandpass filtering (Hz)."},
-            {"name": "PSP_HIGH_FREQ", "default": 1.5, "type": float, "depends_on": "PSP", "advanced": True, "help": "Upper frequency cutoff for PSP bandpass filtering (Hz)."},
+            {"name": "PSP_USE_HEART_RATE_BAND", "default": True, "type": bool, "depends_on": [{"parent_name": "PSP"}, {"parent_name": "HEART_RATE"}], "advanced": False, "help": "Adjust the PSP frequency band using the participant's estimated heart rate."},
+            {"name": "PSP_LOW_FREQ", "default": 0.7, "type": float, "depends_on": "PSP", "advanced": True, "help": "Lower frequency bound for the signal band used in PSP calculation (Hz)."},
+            {"name": "PSP_HIGH_FREQ", "default": 1.5, "type": float, "depends_on": "PSP", "advanced": True, "help": "Upper frequency bound for the signal band used in PSP calculation (Hz)."},
+            {"name": "PSP_TIME_WINDOW", "default": 3, "type": int, "depends_on": "PSP", "advanced": False, "help": "Duration of each independent PSP calculation window in seconds."},
+            {"name": "PSP_THRESHOLD", "default": 0.1, "type": float, "depends_on": "PSP", "advanced": False, "help": "PSP threshold on a scale of 0-1. Channels below this value will be marked bad."},
         ]
     },
     {