implement the new preferences in the backend

flares.py

@@ -135,27 +135,40 @@ OPTODE_PLACEMENT: bool
 SHOW_OPTODE_NAMES: bool
 
 SHORT_CHANNELS: bool
+LONG_CHANNELS: bool
 SHORT_CHANNELS_THRESHOLD: float
 LONG_CHANNELS_THRESHOLD: float
 
 HEART_RATE: bool
 SECONDS_TO_STRIP_HR: int
+HR_LOW_FREQ: int
+HR_HIGH_FREQ: int
+HR_SEARCH_MIN: int
+HR_SEARCH_MAX: int
 MAX_LOW_HR: int
 MAX_HIGH_HR: int
 SMOOTHING_WINDOW_HR: int
 HEART_RATE_WINDOW: int
 
 SCI: bool
+SCI_USE_HEART_RATE_BAND: bool
+SCI_LOW_FREQ: float
+SCI_HIGH_FREQ: float
 SCI_TIME_WINDOW: int
 SCI_THRESHOLD: float
 
 SNR: bool
-# SNR_TIME_WINDOW : int #TODO: is this needed?
 SNR_THRESHOLD: float
+SNR_SIGNAL_LOW_FREQ: float
+SNR_SIGNAL_HIGH_FREQ: float
+SNR_NOISE_LOW_FREQ: float
+SNR_NOISE_HIGH_FREQ: float
 
 PSP: bool
 PSP_TIME_WINDOW: int
 PSP_THRESHOLD: float
+PSP_LOW_FREQ: float
+PSP_HIGH_FREQ: float
 
 COEFF_VAR: bool
 COEFF_VAR_THRESHOLD: int
@@ -166,6 +179,8 @@ MAD_THRESHOLD: int
 PSD_NOISE: bool
 TARGET_FREQ_DIV: int
 DB_LIMIT: int
+PSD_MIN_FREQ: float
+PSD_TARGET_BANDWIDTH: float
 
 SENSOR_DROPOUT: bool
 SENSOR_DROPOUT_VARIANCE_THRESHOLD: float
@@ -189,27 +204,48 @@ PPF_UPPER_WAVELENGTH: float
 ENHANCE_NEGATIVE_CORRELATION: bool
 
 FILTER: bool
+FILTER_ALGORITHM: list
 L_FREQ: float
 H_FREQ: float
 L_TRANS_BANDWIDTH: float
 H_TRANS_BANDWIDTH: float
+IIR_TYPE: list
+IIR_ORDER: int
+FILTER_LENGTH: str
+FILTER_PHASE: list
+FIR_WINDOW: list
+FIR_DESIGN: list
+IIR_OUTPUT: list
+PASSBAND_RIPPLE: float
+STOPBAND_ATTENUATION: float
+FILTER_PAD: list
+SKIP_BY_ANNOTATION: list
+FILTER_N_JOBS: int
 
+EVENTS: bool
+EVENT_ID: str
+EVENT_REGEX: str
+EVENT_CHUNK_DURATION: float
+
+EPOCHS: bool
 EPOCH_HANDLING: str
 MAX_SHIFT: int
 T_MIN: int
 T_MAX: int
+BASELINE: list
+REJECT_EPOCHS: bool
+REJECT_HBO_THRESHOLD: float
 
 RESAMPLE: bool
 RESAMPLE_FREQ: int
-STIM_DUR: float
 HRF_MODEL: str
+STIM_DUR: float
+FIR_DELAYS: range
 DRIFT_MODEL: str
 HIGH_PASS: float
 DRIFT_ORDER: int
-FIR_DELAYS: range
 MIN_ONSET: int
 OVERSAMPLING: int
-REMOVE_EVENTS: list
 SHORT_CHANNEL_REGRESSION: bool
 
 NOISE_MODEL: str
@@ -828,7 +864,7 @@ def scalp_coupling_index_windowed_raw(data, time_window: float = 3.0, l_freq: fl
 
     return data, scores, times
 
-def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5):
+def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5, time_window: int = 3, threshold: float = 0.6):
     """
     Calculate the scalp coupling index (SCI) and identify bad channels based on a threshold.
 
@@ -852,21 +888,21 @@ def calculate_scalp_coupling(data, l_freq: float = 0.7, h_freq: float = 1.5):
     print("Calculating scalp coupling index...")
 
     # Compute the SCI
-    _, scores, times = scalp_coupling_index_windowed_raw(data, time_window=SCI_TIME_WINDOW, l_freq=l_freq, h_freq=h_freq)
+    _, scores, times = scalp_coupling_index_windowed_raw(data, time_window=time_window, l_freq=l_freq, h_freq=h_freq)
 
     # Identify channels that don't meet the provided threshold
     print("Identifying channels that do not meet the threshold...")
     sci = scores.mean(axis=1)
-    data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), sci < SCI_THRESHOLD))
+    data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), sci < threshold))
 
     # Determine the colors based on the threshold, and create the figures
     print("Creating the figures...")
-    color_stops = ([0.0, SCI_THRESHOLD, SCI_THRESHOLD+0.1, 0.8, 1.0], [0.0, SCI_THRESHOLD, SCI_THRESHOLD, 1.0])
+    color_stops = ([0.0, threshold, threshold+0.1, 0.8, 1.0], [0.0, threshold, threshold, 1.0])
-    fig1, fig2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, SCI_THRESHOLD, "Scalp Coupling Index")
+    fig1, fig2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, threshold, "Scalp Coupling Index")
 
     print("Successfully calculated scalp coupling index.")
 
-    return list(compress(cast(list[str], getattr(data, "ch_names")), sci < SCI_THRESHOLD)), fig1, fig2
+    return list(compress(cast(list[str], getattr(data, "ch_names")), sci < threshold)), fig1, fig2
 
 
 
@@ -933,7 +969,7 @@ def get_hbo_hbr_picks(raw):
     return hbo_picks, hbr_picks, hbo_wl, hbr_wl
 
 
-def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2):
+def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2, short_channels_threshold=0.015):
     """
     Interpolate bad fNIRS channels using a distance-weighted average of nearby good channels.
 
@@ -984,7 +1020,7 @@ def interpolate_fNIRS_bads_weighted_average(raw, max_dist=0.03, min_neighbors=2)
         hbo_ch = raw.ch_names[hbo_picks[i]]
         hbr_ch = raw.ch_names[hbr_picks[i]]
         is_bad = (hbo_ch in raw.info['bads']) or (hbr_ch in raw.info['bads'])
-        is_short = pair_distances[i] < SHORT_CHANNELS_THRESHOLD
+        is_short = pair_distances[i] < short_channels_threshold
 
         if is_bad:
             bad_pairs.append(i)
@@ -1168,7 +1204,7 @@ def calculate_signal_noise_ratio(data):
 
 
 
-def calculate_peak_power(data: BaseRaw, l_freq: float = 0.7, h_freq: float = 1.5) -> tuple[list[str], Figure, Figure]:
+def calculate_peak_power(data: BaseRaw, time_window: int = 3, threshold: float = 0.1, l_freq: float = 0.7, h_freq: float = 1.5) -> tuple[list[str], Figure, Figure]:
     """
     Calculate peak spectral power (PSP) for fNIRS channels and identify bad channels.
 
@@ -1191,18 +1227,18 @@ def calculate_peak_power(data: BaseRaw, l_freq: float = 0.7, h_freq: float = 1.5
 
     
     # Compute the PSP
-    _, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power(data, time_window=PSP_TIME_WINDOW, threshold=PSP_THRESHOLD, l_freq=l_freq, h_freq=h_freq))
+    _, scores, times = cast(tuple[NDArray[float64], NDArray[float64], list[tuple[float]]], peak_power(data, time_window=time_window, threshold=threshold, l_freq=l_freq, h_freq=h_freq))
 
     # 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 < PSP_THRESHOLD))
+    data.info["bads"] = list(compress(cast(list[str], getattr(data, "ch_names")), psp < threshold))
 
     # Determine the colors based on the threshold, and create the figures
-    color_stops = ([0.0, PSP_THRESHOLD, PSP_THRESHOLD+0.1, PSP_THRESHOLD+0.2, 1.0], [0.0, PSP_THRESHOLD, PSP_THRESHOLD, 1.0])
+    color_stops = ([0.0, threshold, threshold+0.1, threshold+0.2, 1.0], [0.0, threshold, threshold, 1.0])
-    psp1, psp2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, PSP_THRESHOLD, "Peak Spectral Power")
+    psp1, psp2 = plot_timechannel_quality_metrics(data, scores, times, color_stops, threshold, "Peak Spectral Power")
 
 
-    return list(compress(cast(list[str], getattr(data, "ch_names")), psp < PSP_THRESHOLD)), psp1, psp2
+    return list(compress(cast(list[str], getattr(data, "ch_names")), psp < threshold)), psp1, psp2
 
 
 def mark_bads(raw, bad_sci, bad_snr, bad_psp,  bad_coeff_var, bad_range, bad_noise, bad_disp):
@@ -1279,18 +1315,49 @@ def mark_bads(raw, bad_sci, bad_snr, bad_psp,  bad_coeff_var, bad_range, bad_noi
 
 
 
-def filter_the_data(raw_haemo):
+def filter_the_data(
+        raw_haemo,
+        filter_algorithm,
+        l_freq,
+        h_freq,
+        l_trans_bandwidth,
+        h_trans_bandwidth,
+        iir_type,
+        iir_order,
+        filter_length,
+        filter_phase,
+        fir_window,
+        fir_design,
+        iir_output,
+        passband_ripple,
+        stopband_attenuation,
+        filter_pad,
+        skip_by_annotation,
+        filter_n_jobs,
+        verbosity
+    ):
         # --- STEP 5: Filtering (0.01-0.2 Hz bandpass) ---
     fig_filter = raw_haemo.compute_psd(fmax=3).plot(
         average=True, color="r", show=False, amplitude=True
     )
     
-    if L_FREQ == 0 and H_FREQ != 0:
+    if l_freq == 0 and h_freq != 0:
-        raw_haemo = raw_haemo.filter(l_freq=None, h_freq=H_FREQ, h_trans_bandwidth=H_TRANS_BANDWIDTH)
+        if filter_algorithm == "fir":
-    elif L_FREQ != 0 and H_FREQ == 0:
+            raw_haemo = raw_haemo.filter(l_freq=None, h_freq=h_freq, filter_length=filter_length, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, phase=filter_phase, fir_window=fir_window, fir_design=fir_design, pad=filter_pad, verbose=verbosity)
-        raw_haemo = raw_haemo.filter(l_freq=L_FREQ, h_freq=None, l_trans_bandwidth=L_TRANS_BANDWIDTH)
+        else:
-    elif L_FREQ != 0 and H_FREQ != 0:
+            raw_haemo = raw_haemo.filter(l_freq=None, h_freq=h_freq, filter_length=filter_length, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity)
-        raw_haemo = raw_haemo.filter(l_freq=L_FREQ, h_freq=H_FREQ, l_trans_bandwidth=L_TRANS_BANDWIDTH, h_trans_bandwidth=H_TRANS_BANDWIDTH)
+
+    elif l_freq != 0 and h_freq == 0:
+        if filter_algorithm == "fir":
+            raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=None, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, phase=filter_phase, fir_window=fir_window, fir_design=fir_design, pad=filter_pad, verbose=verbosity)
+        else:
+            raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=None, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity)
+
+    elif l_freq != 0 and h_freq != 0:
+        if filter_algorithm == "fir":
+            raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, phase=filter_phase, fir_window=fir_window, fir_design=fir_design, pad=filter_pad, verbose=verbosity)
+        else:
+            raw_haemo = raw_haemo.filter(l_freq=l_freq, h_freq=h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, n_jobs=filter_n_jobs, method=filter_algorithm, pad=filter_pad, verbose=verbosity)
     else:
         print("No filter")
     #raw_haemo = raw_haemo.filter(l_freq=None, h_freq=0.4, h_trans_bandwidth=0.2)
@@ -1307,23 +1374,30 @@ def filter_the_data(raw_haemo):
 
 
 
-def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline):
+def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline, max_shift, reject_epochs, reject_hbo_threshold):
     """
     Attempts to create epochs, shifting event times slightly if 
     sample collisions are detected.
     """
     shift_increment = 1.0 / raw.info['sfreq']  # The duration of exactly one sample
-    #TODO: User expose this
+
-    reject_criteria = dict(hbo=80e-7)
+    for attempt in range(max_shift):  # Limit attempts to avoid infinite loops
-    for attempt in range(MAX_SHIFT):  # Limit attempts to avoid infinite loops
         try:
-            epochs = Epochs(
+            if reject_epochs:
-                raw, events, event_id=event_dict,
+                epochs = Epochs(
-                tmin=tmin, tmax=tmax, baseline=baseline,
+                    raw, events, event_id=event_dict,
-                reject=reject_criteria,
+                    tmin=tmin, tmax=tmax, baseline=baseline,
-                preload=True, verbose=False
+                    reject=reject_hbo_threshold,
-            )
+                    preload=True, verbose=False
+                )
+            else:
+                epochs = Epochs(
+                    raw, events, event_id=event_dict,
+                    tmin=tmin, tmax=tmax, baseline=baseline,
+                    preload=True, verbose=False
+                )
             return epochs
+        
         except RuntimeError as e:
             if "Event time samples were not unique" in str(e):
                 # Find duplicates in the events array (column 0 is the sample index)
@@ -1344,13 +1418,16 @@ def safe_create_epochs(raw, events, event_dict, tmin, tmax, baseline):
 
 
 
-def epochs_calculations(raw_haemo, events, event_dict):
+def epochs_calculations(raw_haemo, events, event_dict, epoch_handling, max_shift, t_min, t_max, baseline, reject_epochs, reject_hbo_threshold):
     fig_epochs = []  # List to store figures
 
-    if EPOCH_HANDLING == 'shift':
+    if epoch_handling == 'shift':
-        epochs = safe_create_epochs(raw=raw_haemo, events=events, event_dict=event_dict, tmin=T_MIN, tmax=T_MAX, baseline=(None, 0))
+        epochs = safe_create_epochs(raw=raw_haemo, events=events, event_dict=event_dict, tmin=t_min, tmax=t_max, baseline=baseline, max_shift=max_shift, reject_epochs=reject_epochs, reject_hbo_threshold=reject_hbo_threshold)
     else:
-        epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=T_MIN, tmax=T_MAX, baseline=(None, 0))
+        if reject_epochs:
+            epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=t_min, tmax=t_max, baseline=baseline, reject=reject_hbo_threshold)
+        else:
+            epochs = Epochs(raw_haemo, events, event_id=event_dict, tmin=t_min, tmax=t_max, baseline=baseline)
 
     # Make a copy of the epochs and drop bad ones
     epochs2 = epochs.copy()
@@ -1502,7 +1579,25 @@ def epochs_calculations(raw_haemo, events, event_dict):
 
 
 
-def make_design_matrix(raw_haemo):
+def make_design_matrix(
+        raw_haemo,
+        resample,
+        resample_freq,
+        stim_dur,
+        hrf_model,
+        drift_model,
+        high_pass,
+        drift_order,
+        fir_delays,
+        min_onset,
+        oversampling,
+        short_channel_regression,
+        short_channels,
+        long_channels,
+        short_channels_threshold,
+        long_channels_threshold,
+        folding_bypass
+        ):
 
     # events_to_remove = REMOVE_EVENTS
     events_to_remove = ""
@@ -1515,37 +1610,39 @@ def make_design_matrix(raw_haemo):
         description=[ann['description'] for ann in filtered_annotations]
     )
     
-    if SHORT_CHANNELS:
+    if short_channels:
-        short_chans = get_short_channels(raw_haemo, max_dist=SHORT_CHANNELS_THRESHOLD)
+        short_chans = get_short_channels(raw_haemo, max_dist=short_channels_threshold)
-        raw_haemo = get_long_channels(raw_haemo, min_dist=SHORT_CHANNELS_THRESHOLD, max_dist=LONG_CHANNELS_THRESHOLD)
+        if long_channels:
+            raw_haemo = get_long_channels(raw_haemo, min_dist=long_channels_threshold, max_dist=long_channels_threshold)
+
     else:
         short_chans = None
 
     # Set the new annotations
     raw_haemo.set_annotations(new_annot)
 
-    if RESAMPLE:
+    if resample:
-        raw_haemo.resample(RESAMPLE_FREQ, npad="auto")
+        raw_haemo.resample(resample_freq, npad="auto")
         raw_haemo._data = raw_haemo._data * 1e6
         try:
-            short_chans.resample(RESAMPLE_FREQ)
+            short_chans.resample(resample_freq)
         except:
             pass
 
     design_matrix = make_first_level_design_matrix(
         raw=raw_haemo,
-        stim_dur=STIM_DUR,
+        stim_dur=stim_dur,
-        hrf_model=HRF_MODEL,
+        hrf_model=hrf_model,
-        drift_model=DRIFT_MODEL,
+        drift_model=drift_model,
-        high_pass=HIGH_PASS,
+        high_pass=high_pass,
-        drift_order=DRIFT_ORDER,
+        drift_order=drift_order,
-        fir_delays=FIR_DELAYS,
+        fir_delays=fir_delays,
-        min_onset=MIN_ONSET,
+        min_onset=min_onset,
-        oversampling=OVERSAMPLING
+        oversampling=oversampling
     )
 
     # 3) Average and Append Short Channels
-    if SHORT_CHANNEL_REGRESSION and not FOLDING_BYP:
+    if short_channel_regression and not folding_bypass:
         if short_chans is not None and len(short_chans.ch_names) > 0:
             ch_types = short_chans.get_channel_types()
             
@@ -2724,7 +2821,7 @@ def plot_fir_model_results(
 
 
 
-def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]:
+def load_snirf(file_path: str, downsample_frequency: int, verbosity: bool) -> tuple[BaseRaw, Figure]:
     """
     Loads a snirf file, optionally drops channels, downsamples, and creates a figure showing the results.
 
@@ -2745,7 +2842,7 @@ def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]:
     """
     
     # Read the snirf file
-    raw = read_raw_snirf(file_path, preload=True, verbose=VERBOSITY) # type: ignore
+    raw = read_raw_snirf(file_path, preload=True, verbose=verbosity) # type: ignore
     #raw.load_data(verbose=VERBOSITY) # type: ignore  redundant since preload is set to true
 
     # TODO: Why was this commented again?
@@ -2772,7 +2869,7 @@ def load_snirf(file_path: str) -> tuple[BaseRaw, Figure]:
     if DOWNSAMPLE:
         logger.info("Downsample was specified.")
         sfreq_old = getattr(raw, "info")["sfreq"]
-        raw.resample(DOWNSAMPLE_FREQUENCY, verbose=VERBOSITY) # type: ignore
+        raw.resample(downsample_frequency, verbose=verbosity) # type: ignore
         sfreq_new = getattr(raw, "info")["sfreq"]
         logger.info(f"Finished downsampling. Old frequency: {sfreq_old}. New frequency: {sfreq_new}.")
 
@@ -4437,7 +4534,7 @@ def iqr_threshold(coeffs: NDArray[float64], k: float = 1.5) -> floating[Any]:
 
 
 
-def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: int = 3) -> NDArray[float64]:
+def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: int = 3, iqr: float = 1.5) -> NDArray[float64]:
     """
     Denoises a signal using wavelet decomposition and IQR-based thresholding on detail coefficients.
 
@@ -4463,7 +4560,7 @@ def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: i
 
     # Threshold detail coefficients to reduce noise
     for cD in coeffs[1:]:
-        threshold = iqr_threshold(cD, IQR)
+        threshold = iqr_threshold(cD, iqr)
         cD_thresh = np.sign(cD) * np.maximum(np.abs(cD) - threshold, 0.0) # np.where((cD < lower) | (cD > upper), 0, cD)
         cD_thresh = cD_thresh.astype(float64)
         denoised_coeffs.append(cD_thresh)
@@ -4474,7 +4571,7 @@ def wavelet_iqr_denoise(signal: NDArray[float64], wavelet: str = 'db4', level: i
         
 
 
-def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]:
+def calculate_and_apply_wavelet(data: BaseRaw, wavelet_type: str, wavelet_level: int, iqr: float, verbosity: bool) -> tuple[BaseRaw, Figure]:
     """
     Applies a wavelet IQR denoising filter to the data and generates a plot.
 
@@ -4496,17 +4593,17 @@ def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]:
 
     # Denoise the data
     logger.info("Denoising the data...")
-    loaded_data: NDArray[float64] = data.get_data(verbose=VERBOSITY) # type: ignore
+    loaded_data: NDArray[float64] = data.get_data(verbose=verbosity) # type: ignore
     denoised_data = np.zeros_like(loaded_data)
 
     logger.info("Calculating the IQR, decomposing the signal, and thresholding the coefficients...")
     for ch in range(loaded_data.shape[0]):
-        denoised_data[ch, :] = wavelet_iqr_denoise(loaded_data[ch, :], wavelet=WAVELET_TYPE, level=WAVELET_LEVEL)
+        denoised_data[ch, :] = wavelet_iqr_denoise(loaded_data[ch, :], wavelet=wavelet_type, level=wavelet_level, iqr=iqr)
 
     # Reconstruct the data with the annotations
     logger.info("Reconstructing the data with annotations...")
-    raw_with_tddr_and_wavelet = RawArray(denoised_data, cast(Info, data.info), verbose=VERBOSITY)
+    raw_with_tddr_and_wavelet = RawArray(denoised_data, cast(Info, data.info), verbose=verbosity)
-    raw_with_tddr_and_wavelet.set_annotations(data.annotations.copy(), verbose=VERBOSITY) # type: ignore
+    raw_with_tddr_and_wavelet.set_annotations(data.annotations.copy(), verbose=verbosity) # type: ignore
 
     # Create a figure for the results
     logger.info("Creating the figure...")
@@ -4521,7 +4618,7 @@ def calculate_and_apply_wavelet(data: BaseRaw) -> tuple[BaseRaw, Figure]:
 
 
 
-def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None) -> tuple[float, NDArray[float64], NDArray[float64]]:
+def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None, seconds_to_strip_hr: int, verbosity: bool) -> tuple[float, NDArray[float64], NDArray[float64]]:
     """
     Extract and trim short-channel fNIRS signal for heart rate analysis.
 
@@ -4546,16 +4643,16 @@ def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None)
     # If a short channel exists, use it for our signal. Otherwise just take the first channel in the data
     # TODO: Find a better way around this
     if short_chans is not None:
-        signal = cast(NDArray[float64], short_chans.get_data(picks=[0], verbose=VERBOSITY))[0] # type: ignore
+        signal = cast(NDArray[float64], short_chans.get_data(picks=[0], verbose=verbosity))[0] # type: ignore
     else:
-        signal = cast(NDArray[float64], data.get_data(picks=[0], verbose=VERBOSITY))[0] # type: ignore
+        signal = cast(NDArray[float64], data.get_data(picks=[0], verbose=verbosity))[0] # type: ignore
 
     # Calculate the sampling frequency
     sfreq = cast(int, data.info['sfreq'])
 
     # Trim start and end of the signal to remove edge artifacts
-    logger.info(f"Removing {SECONDS_TO_STRIP_HR} seconds from the beginning and end of the file...")
+    logger.info(f"Removing {seconds_to_strip_hr} seconds from the beginning and end of the file...")
-    strip_samples = int(sfreq * SECONDS_TO_STRIP_HR)
+    strip_samples = int(sfreq * seconds_to_strip_hr)
     signal_trimmed = signal[strip_samples:-strip_samples]
     times_trimmed = cast(NDArray[float64], getattr(data, "times"))[strip_samples:-strip_samples]
 
@@ -4563,7 +4660,7 @@ def short_channel_processing_for_hr(data: BaseRaw, short_chans: BaseRaw | None)
 
 
 
-def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64]) -> tuple[NDArray[float64], float]:
+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) -> tuple[NDArray[float64], float]:
     """
     Calculate and smooth heart rate from a trimmed signal using NeuroKit.
 
@@ -4585,19 +4682,19 @@ def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64]
 
     # Filter signal to isolate heart rate frequencies and detect peaks
     logger.info("Filtering the signal and detecting peaks...")
-    signal_filtered = cast(NDArray[float64], nk.signal_filter(signal_trimmed, sampling_rate=sfreq, lowcut=0.8, highcut=2.5)) # type: ignore
+    signal_filtered = cast(NDArray[float64], nk.signal_filter(signal_trimmed, sampling_rate=sfreq, lowcut=hr_low_freq, highcut=hr_high_freq)) # type: ignore
     peaks_dict = cast(dict[str, Any], nk.signal_findpeaks(signal_filtered)) # type: ignore
     peaks = peaks_dict['Peaks']
     hr = cast(NDArray[float64], nk.signal_rate(peaks, sampling_rate=sfreq, desired_length=len(signal_trimmed))) # type: ignore
-    hr_clean = np.clip(hr, MAX_LOW_HR, MAX_HIGH_HR)
+    hr_clean = np.clip(hr, max_low_hr, max_high_hr)
 
     # Smooth heart rate time series by replacing spikes with local rolling mean and calculate the mean
     logger.info("Smoothing the signal and calculating the mean...")
     hr_series = pd.Series(hr_clean)
-    local_median = hr_series.rolling(window=SMOOTHING_WINDOW_HR, center=True, min_periods=1).median()
+    local_median = hr_series.rolling(window=smoothing_window_hr, center=True, min_periods=1).median()
     spikes = hr_series > (local_median + 10)
     smoothed_values = hr_series.copy()
-    smoothed_spikes = hr_series.rolling(window=SMOOTHING_WINDOW_HR, center=True, min_periods=1).mean()
+    smoothed_spikes = hr_series.rolling(window=smoothing_window_hr, center=True, min_periods=1).mean()
     smoothed_values[spikes] = smoothed_spikes[spikes]
     hr_smooth_nk = cast(NDArray[float64], smoothed_values.to_numpy()) # type: ignore
     mean_hr_nk = hr_smooth_nk.mean()
@@ -4612,7 +4709,7 @@ def calculate_heart_rate_neurokit(sfreq: float, signal_trimmed: NDArray[float64]
 
 
        
-def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) -> tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]:
+def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64], search_min, search_max) -> tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]:
     """
     Estimate heart rate using spectral analysis on a high-pass filtered signal.
 
@@ -4628,7 +4725,7 @@ def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) -
     tuple[NDArray[floating[Any]], NDArray[float64], np.ndarray[Any, np.dtype[np.bool_]], float]
         - NDArray[floating[Any]]: Frequencies converted to beats per minute (BPM).
         - NDArray[float64]: Power spectral density (PSD) of the signal.
-        - np.ndarray[Any, np.dtype[np.bool_]]: Boolean mask indicating frequencies within heart rate range (30-300 BPM).
+        - np.ndarray[Any, np.dtype[np.bool_]]: Boolean mask indicating frequencies within heart rate range.
         - float: Estimated mean heart rate in BPM corresponding to the PSD peak within the range.
     """
 
@@ -4644,10 +4741,10 @@ def calculate_heart_rate_scipy(sfreq: float, signal_trimmed: NDArray[float64]) -
     nperseg = min(len(signal_hp), 4096)
     frequencies_scipy, psd_scipy = cast(tuple[NDArray[float64], NDArray[float64]], welch(signal_hp, fs=sfreq, nperseg=nperseg, noverlap=nperseg//2))
 
-    # Convert frequency values to beats per minute (BPM) and set a heart rate range (30-300 BPM)
+    # Convert frequency values to beats per minute (BPM) and set a heart rate range
     logger.info("Converting to BPM...")
     freq_bpm_scipy = frequencies_scipy * 60
-    freq_range_scipy = (freq_bpm_scipy > 30) & (freq_bpm_scipy < 300)
+    freq_range_scipy = (freq_bpm_scipy > search_min) & (freq_bpm_scipy < search_max)
 
     # Identify the peak frequency within the heart rate range and estimate the mean heart rate in BPM
     logger.info("Finding the mean...")
@@ -4667,7 +4764,8 @@ def plot_heart_rate(
     hr_smooth_nk: NDArray[floating[Any]],
     mean_hr_nk: float,
     times_trimmed: NDArray[floating[Any]],
-    overruled: bool
+    overruled: bool,
+    hr_window: int
 ) -> tuple[Figure, Figure]:
     """
     Generate plots comparing heart rate estimates from SciPy PSD and NeuroKit2.
@@ -4704,7 +4802,7 @@ def plot_heart_rate(
     ax1.set_xlim(30, 300)
     ax1.plot(freq_bpm_scipy[freq_range_scipy], psd_scipy[freq_range_scipy]) # type: ignore
     ax1.axvline(x=mean_hr_scipy, color='red', linestyle='--', label=f'Mean HR: {mean_hr_scipy:.1f} BPM') # type: ignore
-    ax1.axvspan(min(mean_hr_nk - HEART_RATE_WINDOW, mean_hr_scipy - HEART_RATE_WINDOW), max(mean_hr_nk + HEART_RATE_WINDOW, mean_hr_scipy + HEART_RATE_WINDOW), color='yellow', alpha=0.3, label=f'HR Range ±{HEART_RATE_WINDOW} BPM') # type: ignore
+    ax1.axvspan(min(mean_hr_nk - hr_window, mean_hr_scipy - hr_window), max(mean_hr_nk + hr_window, mean_hr_scipy + hr_window), color='yellow', alpha=0.3, label=f'HR Range ±{hr_window} BPM') # type: ignore
     ax1.set_xlabel('Heart Rate (BPM)') # type: ignore
     ax1.set_ylabel('Power Spectral Density') # type: ignore
     ax1.set_title('PSD of fNIRS signal - Peak indicates Heart Rate') # type: ignore
@@ -4962,7 +5060,7 @@ def detect_sensor_dropout(raw, threshold_ratio=0.05):
 
 
 
-def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4):
+def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4, min_freq=0.1, target_bandwith=0.2):
     """
     Identifies channels with excessive power at high frequencies 
     (sfreq/4), usually indicating electronic interference.
@@ -4972,11 +5070,11 @@ def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4):
     target_freq = sfreq / freq_div
     
     # Compute PSD
-    spectrum = raw.compute_psd(fmin=0.1, fmax=sfreq/2)
+    spectrum = raw.compute_psd(fmin=min_freq, fmax=sfreq/2)
     psd_data, freqs = spectrum.get_data(return_freqs=True)
     
     # Find power near the target frequency
-    f_idx = np.where((freqs >= target_freq - 0.2) & (freqs <= target_freq + 0.2))[0]
+    f_idx = np.where((freqs >= target_freq - target_bandwith) & (freqs <= target_freq + target_bandwith))[0]
     power_at_target = np.mean(psd_data[:, f_idx], axis=1)
     
     abs_threshold = 10 ** (db_limit / 10)
@@ -5003,7 +5101,7 @@ def detect_spectral_noise_spike(raw, db_limit=-60, freq_div=4):
 
 
 
-def find_bad_channels_by_amplitude_range(raw, threshold=4.0):
+def find_bad_channels_by_amplitude_range(raw, threshold=4):
     """Median absolute deviation"""
     picks = [ch for ch in raw.ch_names]
     data = raw.get_data(picks=picks)
@@ -5094,14 +5192,14 @@ def find_bad_channels_coeff_var(raw, coeff_var_threshold=25.0):
 
 
 
-def hr_calc(raw):
+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):
-    if SHORT_CHANNELS:
+    if short_channels:
-        short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD)
+        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)
+    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_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)
-    freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy = calculate_heart_rate_scipy(sfreq, signal_trimmed)
+    freq_bpm_scipy, psd_scipy, freq_range_scipy, mean_hr_scipy = calculate_heart_rate_scipy(sfreq, signal_trimmed, search_min=search_min, search_max=search_max)
     
     # HACK: This sucks but looking at the graphs I trust neurokit2 more
     overruled = False
@@ -5112,7 +5210,7 @@ def hr_calc(raw):
         mean_hr_scipy = mean_hr_nk
         overruled = True
 
-    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)
+    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)
    
     fig = raw.plot_psd(show=False)
     raw_filtered = raw.copy().filter(0.5, 3, fir_design='firwin')
@@ -5123,7 +5221,7 @@ def hr_calc(raw):
     # --- Parameters for PSD ---
     desired_bin_hz = 0.1
     nperseg = int(sfreq / desired_bin_hz)
-    hr_range = (30, 180) # TODO: Should this not use the user defined values?
+    hr_range = (search_min, search_max)
 
     # --- Function to find strongest local peak ---
     def find_hr_from_psd(ch_data):
@@ -5153,11 +5251,11 @@ def hr_calc(raw):
     return fig, hr1, hr2, low, high
 
 
-def trim_participant_data(raw):
+def trim_participant_data(raw, seconds_to_keep: float):
     if hasattr(raw, 'annotations') and len(raw.annotations) > 0:
         # Get time of first event
         first_event_time = raw.annotations.onset[0]
-        trim_time = max(0, first_event_time - SECONDS_TO_KEEP)  # Ensure we don't go negative
+        trim_time = max(0, first_event_time - seconds_to_keep)  # Ensure we don't go negative
         raw.crop(tmin=trim_time)
         # Shift annotation onsets to match new t=0
 
@@ -5180,14 +5278,14 @@ def trim_participant_data(raw):
     return raw, fig_trimmed
 
 
-def remove_bad_channels(raw, bad_channels):
+def remove_bad_channels(raw, bad_channels, max_bad_channels: int):
     num_bad = len(bad_channels)
     
     # Check against the threshold
-    if num_bad > MAX_BAD_CHANNELS:
+    if num_bad > max_bad_channels:
         raise Exception(
             f"Data Quality Error: {num_bad} channels flagged for removal, "
-            f"which exceeds the limit of {MAX_BAD_CHANNELS}. To avoid this, "
+            f"which exceeds the limit of {max_bad_channels}. To avoid this, "
             f"either lower your filtering parameters or increase MAX_BAD_CHANNELS."
         )
 
@@ -5196,9 +5294,9 @@ def remove_bad_channels(raw, bad_channels):
     return raw
 
 
-def make_and_run_glm(raw_haemo, df_design_matrix):
+def make_and_run_glm(raw_haemo, df_design_matrix, noise_model, bins, n_jobs, verbosity):
 
-    glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbose=VERBOSITY)
+    glm_est = run_glm(raw_haemo, df_design_matrix, noise_model=noise_model, bins=bins, n_jobs=n_jobs, verbose=verbosity)
     fir_cols = [col for col in df_design_matrix.columns if "_delay_" in col]
     
     if fir_cols:
@@ -5265,11 +5363,11 @@ def generate_channel_results(glm_est, file_path):
     return df_cha
 
 
-def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path):
+def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_location):
     
     rois_formatted = {}
     try:
-        with open(JSON_LOCATION, 'r') as f:
+        with open(json_location, 'r') as f:
             roi_data = json.load(f)
     
         for region in roi_data.get("regions_of_interest", []):
@@ -5305,7 +5403,7 @@ def generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path):
     # ---------------------------------------------------------------------
     if not rois_formatted:
         logger.error(
-            f"'{JSON_LOCATION}' produced zero valid ROIs — attempting automatic "
+            f"'{json_location}' produced zero valid ROIs — attempting automatic "
             f"geometric fallback (Left/Right, then Front/Back)."
         )
         rois_formatted = _build_geometric_fallback_rois(raw_haemo)
@@ -5449,6 +5547,7 @@ def haemoglobin_concentration(raw_od, file_path, override_ppf=False, ppf_lower_w
 
 
 def process_participant(file_path, progress_callback=None):
+    """Goal is to have all of the constants be located here, and nowhere else in the code to remove any ambiguity."""
 
     # Step 0: Setting up
     fig_individual: dict[str, Figure] = {}
@@ -5459,7 +5558,7 @@ def process_participant(file_path, progress_callback=None):
     }
 
     # Step 1: Preprocessing
-    raw = load_snirf(file_path)
+    raw = load_snirf(file_path=file_path, downsample_frequency=DOWNSAMPLE_FREQUENCY, verbosity=VERBOSITY)
     fig_raw = raw.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Loaded Raw", show=False)
     fig_individual["Loaded Raw Data"] = fig_raw
     if progress_callback: progress_callback(1)
@@ -5467,7 +5566,7 @@ def process_participant(file_path, progress_callback=None):
     
     # Step 2: Trimming
     if TRIM and not FOLDING_BYP:
-        raw, fig_trimmed = trim_participant_data(raw)
+        raw, fig_trimmed = trim_participant_data(raw, seconds_to_keep=SECONDS_TO_KEEP)
         fig_individual["Trimmed Raw Data"] = fig_trimmed
     if progress_callback: progress_callback(2)
     logger.info("Step 2 Completed.")
@@ -5485,13 +5584,28 @@ def process_participant(file_path, progress_callback=None):
         _short_chans = get_short_channels(raw, max_dist=SHORT_CHANNELS_THRESHOLD) # JUST FOR PLOTTING THEM SEPERATELY
         fig_short_chans = _short_chans.plot(duration=raw.times[-1], n_channels=raw.info['nchan'], title="Short Channels Only", show=False, verbose=VERBOSITY)
         fig_individual["Short Channels Raw Data"] = fig_short_chans
-    raw = get_long_channels(raw, min_dist=0, max_dist=LONG_CHANNELS_THRESHOLD)
+    if LONG_CHANNELS:
+        raw = get_long_channels(raw, min_dist=0, max_dist=LONG_CHANNELS_THRESHOLD)
     if progress_callback: progress_callback(4)
     logger.info("Step 4 Completed.")
 
     # Step 5: Heart Rate
     if HEART_RATE and not FOLDING_BYP:
-        fig, hr1, hr2, low, high = hr_calc(raw)
+        fig, hr1, hr2, low, high = hr_calc(
+            raw, 
+            seconds_to_strip_hr=SECONDS_TO_STRIP_HR, 
+            l_freq=HR_LOW_FREQ, 
+            h_freq=HR_HIGH_FREQ, 
+            search_min=HR_SEARCH_MIN, 
+            search_max=HR_SEARCH_MAX, 
+            max_low_hr=MAX_LOW_HR, 
+            max_high_hr=MAX_HIGH_HR, 
+            smoothing_window_hr=SMOOTHING_WINDOW_HR, 
+            hr_window=HEART_RATE_WINDOW,
+            short_channels=SHORT_CHANNELS, 
+            short_channels_threshold=SHORT_CHANNELS_THRESHOLD,
+            verbosity=VERBOSITY
+        )
         fig_individual["Power Spectral Density"] = fig
         fig_individual['Heart Rate - PSD'] = hr1
         fig_individual['Heart Rate - Time'] = hr2
@@ -5501,10 +5615,10 @@ def process_participant(file_path, progress_callback=None):
     # Step 6: Scalp Coupling Index    
     bad_sci = []
     if SCI and not FOLDING_BYP:
-        if HEART_RATE:
+        if HEART_RATE and SCI_USE_HEART_RATE_BAND:
-            bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, low, high)
+            bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw, l_freq=low, h_freq=high, time_window=SCI_TIME_WINDOW, threshold=SCI_THRESHOLD)
         else:
-            bad_sci, fig_sci_1, fig_sci_2 = calculate_scalp_coupling(raw)
+            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)
         fig_individual["Scalp Coupling Index Heatmap"] = fig_sci_1
         fig_individual["Scalp Coupling Index Binary Heatmap"] = fig_sci_2
     if progress_callback: progress_callback(6)
@@ -5521,7 +5635,7 @@ def process_participant(file_path, 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)
+        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)
         fig_individual["Peak Spectral Power Heatmap"] = fig_psp1
         fig_individual["Peak Spectral Power Binary Heatmap"] = fig_psp2
     if progress_callback: progress_callback(8)
@@ -5546,7 +5660,7 @@ def process_participant(file_path, progress_callback=None):
     # Step 11: Power Spectral Density Noise
     bad_noise = []
     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)
+        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)
         fig_individual['Power Spectral Density Noise'] = fig_noise
     if progress_callback: progress_callback(11)
     logger.info("Step 11 Completed.")
@@ -5567,11 +5681,11 @@ def process_participant(file_path, progress_callback=None):
             fig_individual["Bad Channels Data"] = fig_raw_before
         if bad_channels:
             if BAD_CHANNELS_HANDLING == "Interpolate":
-                raw, fig_raw_after, fig_compare = interpolate_fNIRS_bads_weighted_average(raw, max_dist=MAX_DIST, min_neighbors=MIN_NEIGHBORS)
+                raw, fig_raw_after, fig_compare = interpolate_fNIRS_bads_weighted_average(raw, max_dist=MAX_DIST, min_neighbors=MIN_NEIGHBORS, short_channels_threshold=SHORT_CHANNELS_THRESHOLD)
                 fig_individual["Data after Interpolating Bad Channels"] = fig_raw_after
                 fig_individual["Bad Channels Interpolation Results"] = fig_compare
             elif BAD_CHANNELS_HANDLING == "Remove":
-                raw = remove_bad_channels(raw, bad_channels)
+                raw = remove_bad_channels(raw, bad_channels, max_bad_channels=MAX_BAD_CHANNELS)
     if progress_callback: progress_callback(13)
     logger.info("Step 13 Completed.")
 
@@ -5592,7 +5706,7 @@ def process_participant(file_path, progress_callback=None):
     
     # Step 16: Wavelet Filtering
     if WAVELET and not FOLDING_BYP:
-        raw_od, fig = calculate_and_apply_wavelet(raw_od)
+        raw_od, fig = calculate_and_apply_wavelet(data=raw_od, wavelet_type=WAVELET_TYPE, wavelet_level=WAVELET_LEVEL, iqr=IQR, verbosity=VERBOSITY)
         fig_individual["Wavelet"] = fig
     if progress_callback: progress_callback(16)
     logger.info("Step 16 Completed.")
@@ -5614,36 +5728,85 @@ def process_participant(file_path, progress_callback=None):
     
     # Step 19: Filter
     if FILTER and not FOLDING_BYP:
-        raw_haemo, fig_filter, fig_raw_haemo_filter = filter_the_data(raw_haemo)
+        raw_haemo, fig_filter, fig_raw_haemo_filter = filter_the_data(
+            raw_haemo,
+            filter_algorithm=FILTER_ALGORITHM,
+            l_freq=L_FREQ,
+            h_freq=H_FREQ,
+            l_trans_bandwidth=L_TRANS_BANDWIDTH,
+            h_trans_bandwidth=H_TRANS_BANDWIDTH,
+            iir_type=IIR_TYPE,
+            iir_order=IIR_ORDER,
+            filter_length=FILTER_LENGTH,
+            filter_phase=FILTER_PHASE,
+            fir_window=FIR_WINDOW,
+            fir_design=FIR_DESIGN,
+            iir_output=IIR_OUTPUT,
+            passband_ripple=PASSBAND_RIPPLE,
+            stopband_attenuation=STOPBAND_ATTENUATION,
+            filter_pad=FILTER_PAD,
+            skip_by_annotation=SKIP_BY_ANNOTATION,
+            filter_n_jobs=FILTER_N_JOBS,
+            verbosity=VERBOSITY
+        )
         fig_individual["Filter_1"] = fig_filter
         fig_individual["Filter_2"] = fig_raw_haemo_filter
     if progress_callback: progress_callback(19)
     logger.info("Step 19 Completed.")
 
     # Step 20: Extracting Events
-    if not FOLDING_BYP:
+    if EVENTS and not FOLDING_BYP:
-        events, event_dict = events_from_annotations(raw_haemo)
+        events, event_dict = events_from_annotations(raw_haemo, event_id=EVENT_ID, regexp=EVENT_REGEX, verbose=VERBOSITY) #TODO: Implement the chunk duration
         fig_events = plot_events(events, event_id=event_dict, sfreq=raw_haemo.info["sfreq"], show=False)
         fig_individual["Events"] = fig_events
     if progress_callback: progress_callback(20)
     logger.info("Step 20 Completed.")
 
     # Step 21: Epoch Calculations
-    if not FOLDING_BYP:
+    if EPOCHS and EVENTS and not FOLDING_BYP:
-        epochs, fig_epochs = epochs_calculations(raw_haemo, events, event_dict)
+        epochs, fig_epochs = epochs_calculations(
+            raw_haemo, 
+            events, 
+            event_dict, 
+            epoch_handling=EPOCH_HANDLING,
+            max_shift=MAX_SHIFT,
+            t_min=T_MIN,
+            t_max=T_MAX,
+            baseline=(None,0), #TODO: Unhardcode this
+            reject_epochs=REJECT_EPOCHS,
+            reject_hbo_threshold=dict(hbo=REJECT_HBO_THRESHOLD)
+        )
         for name, fig in fig_epochs:
             fig_individual[f"epochs_{name}"] = fig
     if progress_callback: progress_callback(21)
     logger.info("Step 21 Completed.")
 
     # Step 22: Design Matrix
-    raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix(raw_haemo) # Short channel is re-applied inside this method
+    raw_haemo, df_design_matrix, fig_design_matrix = make_design_matrix(
+        raw_haemo=raw_haemo,
+        resample=RESAMPLE,
+        resample_freq=RESAMPLE_FREQ,
+        stim_dur=STIM_DUR,
+        hrf_model=HRF_MODEL,
+        drift_model=DRIFT_MODEL,
+        high_pass=HIGH_PASS,
+        drift_order=DRIFT_ORDER,
+        fir_delays=FIR_DELAYS,
+        min_onset=MIN_ONSET,
+        oversampling=OVERSAMPLING,
+        short_channel_regression=SHORT_CHANNEL_REGRESSION,
+        short_channels=SHORT_CHANNELS,
+        long_channels=LONG_CHANNELS,
+        short_channels_threshold=SHORT_CHANNELS_THRESHOLD,
+        long_channels_threshold=LONG_CHANNELS_THRESHOLD,
+        folding_bypass=FOLDING_BYP
+    )
     fig_individual["Design Matrix"] = fig_design_matrix
     if progress_callback: progress_callback(22)
     logger.info("Step 22 Completed.")
 
     # Step 23: General Linear Model
-    glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix)
+    glm_est, fig_glm_topo = make_and_run_glm(raw_haemo, df_design_matrix, noise_model=NOISE_MODEL, bins=BINS, n_jobs=N_JOBS, verbosity=VERBOSITY)
     fig_individual["GLM Topography"] = fig_glm_topo
     if progress_callback: progress_callback(23)
     logger.info("23")
@@ -5662,7 +5825,7 @@ def process_participant(file_path, progress_callback=None):
     logger.info("25")
 
     # Step 26: Generate Region of Interest Results
-    df_roi, fig_roi = generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path)
+    df_roi, fig_roi = generate_roi_results(raw_haemo, df_design_matrix, glm_est, file_path, json_location=JSON_LOCATION)
     fig_individual["Region of Interest"] = fig_roi
     if progress_callback: progress_callback(26)
     logger.info("26")

main.py

@@ -109,7 +109,7 @@ SECTIONS = [
             {"name": "SHORT_CHANNELS", "default": True, "type": bool, "advanced": False, "help": "This should be set to True if the data has a short channel present in the data. For more information about short channels, please visit the Wiki."},
             {"name": "LONG_CHANNELS", "default": True, "type": bool, "advanced": True, "help": "Should channels exceeding the maximum allowed distance be removed?"},
             {"name": "SHORT_CHANNELS_THRESHOLD", "default": 0.015, "type": float, "depends_on": "SHORT_CHANNELS", "advanced": False, "help": "The maximum distance the short channel can be in metres before it is no longer considered a short channel."},
-            {"name": "LONG_CHANNELS_THRESHOLD", "default": 0.045, "type": float, "advanced": False, "help": "The maximum distance channels can be in metres. Any channel longer than this distance will be discarded."},
+            {"name": "LONG_CHANNELS_THRESHOLD", "default": 0.045, "type": float, "depends_on": "LONG_CHANNELS", "advanced": False, "help": "The maximum distance channels can be in metres. Any channel longer than this distance will be discarded."},
         ]
     },
     {
@@ -131,7 +131,7 @@ SECTIONS = [
         "title": "Scalp Coupling Index",
         "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": "SCI", "advanced": False, "help": "Adjust the SCI frequency band using the participant's estimated heart rate."},
+            {"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_TIME_WINDOW", "default": 3, "type": int, "depends_on": "SCI", "advanced": False, "help": "Duration of each independent SCI calculation window in seconds."},
@@ -147,7 +147,6 @@ SECTIONS = [
             {"name": "SNR_SIGNAL_HIGH_FREQ", "default": 0.5, "type": float, "depends_on": "SNR", "advanced": True, "help": "Upper frequency bound for the signal band used in SNR calculation (Hz)."},
             {"name": "SNR_NOISE_LOW_FREQ", "default": 1.0, "type": float, "depends_on": "SNR", "advanced": True, "help": "Lower frequency bound for the noise band used in SNR calculation (Hz)."},
             {"name": "SNR_NOISE_HIGH_FREQ", "default": 10.0, "type": float, "depends_on": "SNR", "advanced": True, "help": "Upper frequency bound for the noise band used in SNR calculation (Hz)."},
-            {"name": "SNR_PLOT_MAX", "default": 20.0, "type": float, "depends_on": "SNR", "advanced": True, "help": "Maximum SNR value displayed on the plot scale (dB)."},
         ]
     },
     {
@@ -216,9 +215,9 @@ SECTIONS = [
         "title": "Wavelet filtering",
         "params": [
             {"name": "WAVELET", "default": True, "type": bool, "advanced": False, "help": "Apply Wavelet filtering. It is a method to filter involving decomposition, threholding, and reconstruction."},
-            {"name": "IQR", "default": 1.5, "type": float, "depends_on": "WAVELET", "advanced": False, "help": "Scaling factor for the Inter-Quartile Range."},
             {"name": "WAVELET_TYPE", "default": "db4", "type": str, "depends_on": "WAVELET", "advanced": False, "help": "Wavelet type. Valid values are ['bior1.1', 'bior1.3', 'bior1.5', 'bior2.2', 'bior2.4', 'bior2.6', 'bior2.8', 'bior3.1', 'bior3.3', 'bior3.5', 'bior3.7', 'bior3.9', 'bior4.4', 'bior5.5', 'bior6.8', 'coif1', 'coif2', 'coif3', 'coif4', 'coif5', 'coif6', 'coif7', 'coif8', 'coif9', 'coif10', 'coif11', 'coif12', 'coif13', 'coif14', 'coif15', 'coif16', 'coif17', 'db1', 'db2', 'db3', 'db4', 'db5', 'db6', 'db7', 'db8', 'db9', 'db10', 'db11', 'db12', 'db13', 'db14', 'db15', 'db16', 'db17', 'db18', 'db19', 'db20', 'db21', 'db22', 'db23', 'db24', 'db25', 'db26', 'db27', 'db28', 'db29', 'db30', 'db31', 'db32', 'db33', 'db34', 'db35', 'db36', 'db37', 'db38', 'dmey', 'haar', 'rbio1.1', 'rbio1.3', 'rbio1.5', 'rbio2.2', 'rbio2.4', 'rbio2.6', 'rbio2.8', 'rbio3.1', 'rbio3.3', 'rbio3.5', 'rbio3.7', 'rbio3.9', 'rbio4.4', 'rbio5.5', 'rbio6.8', 'sym2', 'sym3', 'sym4', 'sym5', 'sym6', 'sym7', 'sym8', 'sym9', 'sym10', 'sym11', 'sym12', 'sym13', 'sym14', 'sym15', 'sym16', 'sym17', 'sym18', 'sym19', 'sym20']"},
             {"name": "WAVELET_LEVEL", "default": 3, "type": int, "depends_on": "WAVELET", "advanced": False, "help": "Wavelet Decomposition level (must be >= 0)."},
+            {"name": "IQR", "default": 1.5, "type": float, "depends_on": "WAVELET", "advanced": False, "help": "Scaling factor for the Inter-Quartile Range."},
         ]
     },
     {
@@ -239,43 +238,43 @@ SECTIONS = [
         "title": "Filtering",
         "params": [
             {"name": "FILTER", "default": True, "type": bool, "advanced": False, "help": "Should the data be bandpass filtered?"},
-            {"name": "FILTER_ALGORITHM", "default": ["FIR"], "type": list, "options": ["FIR", "IIR"], "exclusive": True, "advanced": False, "help": "Filtering algorithm."},
+            {"name": "FILTER_ALGORITHM", "default": ["fir"], "type": list, "options": ["fir", "iir"], "exclusive": True, "advanced": False, "help": "Filtering algorithm."},
             {"name": "L_FREQ", "default": 0.005, "type": float, "depends_on": "FILTER", "advanced": False, "help": "Any frequencies lower than this value will be removed."},
             {"name": "H_FREQ", "default": 0.3, "type": float, "depends_on": "FILTER", "advanced": False, "help": "Any frequencies higher than this value will be removed."},
             {"name": "L_TRANS_BANDWIDTH", "default": 0.002, "type": float, "depends_on": "FILTER", "advanced": True, "help": "Width of the lower transition band to prevent abrupt filter cutoff."},
             {"name": "H_TRANS_BANDWIDTH", "default": 0.002, "type": float, "depends_on": "FILTER", "advanced": True, "help": "Width of the upper transition band to prevent abrupt filter cutoff."},
-            {"name": "IIR_TYPE", "default": ["butterworth"], "type": list, "options": ["butterworth", "chebyshev1", "chebyshev2", "elliptic", "bessel"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "IIR filter design."},
+            # {"name": "IIR_TYPE", "default": ["butterworth"], "type": list, "options": ["butterworth", "chebyshev1", "chebyshev2", "elliptic", "bessel"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "IIR filter design."},
-            {"name": "IIR_ORDER", "default": 4, "type": int, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Order of the IIR filter."},
+            # {"name": "IIR_ORDER", "default": 4, "type": int, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Order of the IIR filter."},
             {"name": "FILTER_LENGTH", "default": "auto", "type": str, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Length of the FIR filter. 'auto' allows automatic selection."},
-            {"name": "FILTER_PHASE", "default": ["zero"], "type": list, "options": ["zero", "zero-double", "minimum", "minimum-half"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Phase response of the FIR filter."},
+            {"name": "FILTER_PHASE", "default": ["zero"], "type": list, "options": ["zero", "zero-double", "minimum", "minimum-half", "linear"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Phase response of the FIR filter."},
             {"name": "FIR_WINDOW", "default": ["hamming"], "type": list, "options": ["hamming", "hann", "blackman"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Window function used when designing the FIR filter."},
             {"name": "FIR_DESIGN", "default": ["firwin"], "type": list, "options": ["firwin", "firwin2"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "FIR", "advanced": True, "help": "Method used to design the FIR filter."},
-            {"name": "IIR_OUTPUT", "default": ["sos"], "type": list, "options": ["sos", "ba", "zpk"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Representation used for IIR filter coefficients."},
+            # {"name": "IIR_OUTPUT", "default": ["sos"], "type": list, "options": ["sos", "ba", "zpk"], "exclusive": True, "depends_on": "FILTER_ALGORITHM", "depends_value": "IIR", "advanced": True, "help": "Representation used for IIR filter coefficients."},
-            {"name": "PASSBAND_RIPPLE", "default": 1.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev1", "elliptic"], "advanced": True, "help": "Maximum allowed ripple in the passband (dB)."},
+            # {"name": "PASSBAND_RIPPLE", "default": 1.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev1", "elliptic"], "advanced": True, "help": "Maximum allowed ripple in the passband (dB)."},
-            {"name": "STOPBAND_ATTENUATION", "default": 40.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev2", "elliptic"], "advanced": True, "help": "Minimum attenuation in the stopband (dB)."},
+            # {"name": "STOPBAND_ATTENUATION", "default": 40.0, "type": float, "depends_on": "IIR_TYPE", "depends_value": ["chebyshev2", "elliptic"], "advanced": True, "help": "Minimum attenuation in the stopband (dB)."},
             {"name": "FILTER_PAD", "default": ["reflect_limited"], "type": list, "options": ["reflect_limited", "reflect", "edge", "constant"], "exclusive": True, "depends_on": "FILTER", "advanced": True, "help": "Padding strategy used during filtering to reduce edge artifacts."},
-            {"name": "SKIP_BY_ANNOTATION", "default": ["edge", "bad_acq_skip"], "type": list, "depends_on": "FILTER", "advanced": True, "help": "Annotations that should be skipped when applying the filter."},
+            # {"name": "SKIP_BY_ANNOTATION", "default": ["edge", "bad_acq_skip"], "type": list, "depends_on": "FILTER", "advanced": True, "help": "Annotations that should be skipped when applying the filter."},
             {"name": "FILTER_N_JOBS", "default": 1, "type": int, "advanced": True, "help": "Number of parallel jobs used during filtering. Use -1 to use all available CPUs."},
         ]
     },
     {
-        "title": "Extracting Events*",
+        "title": "Extracting Events",
         "params": [
-            {"name": "EVENTS", "default": True, "type": bool, "advanced": False, "help": "Extract events from annotations for visualization and downstream event-based analysis."},
+            {"name": "EVENTS", "default": True, "type": bool, "advanced": True, "help": "Extract events from annotations for visualization and downstream event-based analysis."},
             {"name": "EVENT_ID", "default": "auto", "type": str, "advanced": True, "help": "Controls how annotation descriptions are converted into event identifiers. Use 'auto' for automatic event detection."},
-            {"name": "EVENT_REGEX", "default": "^(?![Bb][Aa][Dd]|[Ee][Dd][Gg][Ee]).*$", "type": str, "advanced": True, "help": "Regular expression used to select which annotations are converted into events. By default, bad and edge annotations are ignored."},
+            {"name": "EVENT_REGEX", "default": r"^(?![Bb][Aa][Dd]|[Ee][Dd][Gg][Ee]).*$", "type": str, "advanced": True, "help": "Regular expression used to select which annotations are converted into events. By default, bad and edge annotations are ignored."},
-            {"name": "EVENT_CHUNK_DURATION", "default": None, "type": float, "advanced": True, "help": "If provided, creates repeated events at this interval within longer annotations instead of only using annotation onset times."},
+            # {"name": "EVENT_CHUNK_DURATION", "default": 0.0, "type": float, "advanced": True, "help": "If provided, creates repeated events at this interval within longer annotations instead of only using annotation onset times."},
         ]
     },
         {
         "title": "Epoch Calculations",
         "params": [
-            {"name": "EPOCHS", "default": True, "type": bool, "advanced": False, "help": "Create epochs around extracted events for condition-based analysis."},
+            {"name": "EPOCHS", "default": True, "type": bool, "depends_on": "EVENTS", "advanced": True, "help": "Create epochs around extracted events for condition-based analysis."},
             {"name": "EPOCH_HANDLING", "default": ["shift"], "type": list, "options": ["shift", "strict", "drop"], "exclusive": True, "advanced": False, "help": "How to handle events occurring at the same sample. Shift moves conflicting events forward, strict raises an error, and drop removes conflicting events."},
             {"name": "MAX_SHIFT", "default": 5, "type": int, "depends_on": "EPOCH_HANDLING", "depends_value": "shift", "advanced": True, "help": "Maximum number of samples to shift conflicting events before failing."},
             {"name": "T_MIN", "default": -5.0, "type": float, "advanced": False, "help": "Time in seconds before each event to include in the epoch."},
             {"name": "T_MAX", "default": 15.0, "type": float, "advanced": False, "help": "Time in seconds after each event to include in the epoch."},
-            {"name": "BASELINE", "default": ["pre_event"], "type": list, "options": ["none", "pre_event"], "exclusive": True, "advanced": False, "help": "Baseline correction applied to epochs. Pre-event uses the period before the event as baseline."},
+            # {"name": "BASELINE", "default": ["pre_event"], "type": list, "options": ["none", "pre_event"], "exclusive": True, "advanced": False, "help": "Baseline correction applied to epochs. Pre-event uses the period before the event as baseline."},
             {"name": "REJECT_EPOCHS", "default": True, "type": bool, "advanced": False, "help": "Automatically reject epochs containing excessively large haemoglobin amplitude changes."},
             {"name": "REJECT_HBO_THRESHOLD", "default": 80e-7, "type": float, "depends_on": "REJECT_EPOCHS", "advanced": True, "help": "Maximum allowed HbO amplitude before an epoch is rejected."},
         ]
@@ -310,13 +309,6 @@ SECTIONS = [
             {"name": "JSON_LOCATION", "default": "", "type": "json_file", "advanced": False, "help": "Location of the JSON file containing region of interest results for significance calculations."},
         ]
     },
-    {
-        "title": "Contrast",
-        "params": [
-            {"name": "CONTRAST_BASELINE", "default": True, "type": bool, "advanced": True, "help": "Calculate contrasts comparing each condition against zero/baseline. This may not always represent a meaningful physiological baseline in fNIRS."},
-            {"name": "CONTRAST_PAIRWISE", "default": True, "type": bool, "advanced": True, "help": "Calculate contrasts comparing each experimental condition against every other condition. The number of contrasts increases with the number of conditions."}
-        ]
-    },
     {
         "title": "Other",
         "params": [

src/shared/flaresbasewidget.py

@@ -417,10 +417,20 @@ class ParamSection(QWidget):
                 widget.textChanged.connect(lambda val, p=param["name"]: self.check_if_changed(p, val))
 
             if "depends_on" in param:
+                deps_list = []
+                
+                if isinstance(param["depends_on"], list):
+                    deps_list = param["depends_on"]
+                
+                else:
+                    deps_list.append({
+                        "parent_name": param["depends_on"],
+                        "depends_value": param.get("depends_value", "True")
+                    })
+
                 self.dependencies.append({
                     "child_name": param["name"],
-                    "parent_name": param["depends_on"],
+                    "conditions": deps_list
-                    "depends_value": param.get("depends_value", "True")
                 })
 
             widget.setToolTip(help_text)
@@ -532,23 +542,44 @@ class ParamSection(QWidget):
         """Disables/Enables widgets based on parent selection values."""
         for dep in self.dependencies:
             child_info = self.widgets.get(dep["child_name"])
-            parent_info = self.widgets.get(dep["parent_name"])
+            if not child_info:
+                continue
+
+            # Default to enabled until a condition fails
+            all_conditions_met = True
+
+            for cond in dep["conditions"]:
+                parent_name = cond.get("parent_name") or cond.get("parent")
+                required_value = str(cond.get("depends_value") if "depends_value" in cond else cond.get("value", "True"))
+                
+                parent_info = self.widgets.get(parent_name)
+                if not parent_info:
+                    all_conditions_met = False
+                    break
 
-            if child_info and parent_info:
                 parent_widget = parent_info["widget"]
 
-                # Get current value of parent (works for both bool-combos and list-combos)
+                # Extract current parent value based on widget type
-                current_parent_value = parent_widget.currentText()
+                if isinstance(parent_widget, QComboBox):
+                    current_parent_value = parent_widget.currentText()
+                elif isinstance(parent_widget, QLineEdit):
+                    current_parent_value = parent_widget.text()
+                elif isinstance(parent_widget, QSpinBox):
+                    current_parent_value = str(parent_widget.value())
+                else:
+                    current_parent_value = str(parent_widget)
 
-                # Check if it matches the required value
+                # If any condition fails, flag as false
-                is_active = (current_parent_value == dep["depends_value"])
+                if current_parent_value != required_value:
+                    all_conditions_met = False
+                    break
                 
-                # Toggle the entire row (Button, Label, and Input)
+            # Toggle the entire row (Button, Label, and Input)
-                h_layout = child_info["h_layout"]
+            h_layout = child_info["h_layout"]
-                for i in range(h_layout.count()):
+            for i in range(h_layout.count()):
-                    item = h_layout.itemAt(i).widget()
+                item = h_layout.itemAt(i).widget()
-                    if item:
+                if item:
-                        item.setEnabled(is_active)
+                    item.setEnabled(all_conditions_met)
 
     def _create_multiselect_dropdown(self, items):
         combo = FullClickComboBox()