Add Envelope + Rise Rate SNR estimation metric

[ksseverson57]

Summary
Adds EnvelopeRRSNR, an alternative SNR estimator for preprocessed dF/F fluorescence traces, drop-in compatible with existing derivative-based estimate_snr (see usage below).

Approach
Instead of estimating noise from the sample-to-sample derivative and signal from raw peak heights, this method:

(1) Tracks a valley-anchored tonic envelope beneath the trace, robust to slow drift/bleaching.
(2) Subtracts the tonic envelope to obtain a phasic residual.
(3) Estimates the noise floor robustly from the below-median half of the residual (IQR-based sigma), minimizing effect of putative phasic events.
(4) Detects candidate phasic peaks above a configurable threshold and refines them with a rise-rate gate to reject slow drifts/decay and keep only fast-onset transients.
(5) Reports SNR as a suprathreshold peak-amplitude statistic (median by default or 95th percentile) divided by the estimated noise floor.

Validation
Benchmarked against the derivative-based method across a synthetic SNR sweep (slow baseline shifts, Gaussian noise, skewed fluorescence events).

See benchmark notebook (BWNM_Signal_Quality_Benchmark.ipynb) in this Code Ocean capsule: [https://codeocean.allenneuraldynamics.org/capsule/1266302/tree]

Results:
Detects peak events at least as well as the derivative method, with slightly better accuracy (e.g., precision, recall, and F1 statistic).
Shows a small, consistent linear bias (slope ~1) in its SNR estimate, for which an optional bias correction is provided (bias_correction param / fit_bias_correction_from_benchmark). Env+RR overestimates SNR by a static offset of ~6, so can easily be subtracted to be conservative.
In contrast, the derivative method appears to be highly nonlinear, overshooting at low SNR and saturating at high SNR, with poor noise and signal estimates.

API
EnvelopeRRSNR(fps=20.0, peak_threshold_sd=1.5, signal_statistic='median', bias_correction=None, config=None)
.fit(trace) → EnvelopeRRResult (with snr, noise, peaks, tonic, residual, etc.)
.estimate(trace) → (snr, noise, peaks) tuple, drop-in compatible with the existing estimate_snr(trace, fps) function
.decompose(trace) → (tonic, phasic) component arrays

[hagikent]

looks good to me
I would ask @j-friedrich to have a look and suggest potential improvement

[alexpiet]

Add bwnm_signal_utils, which should be named something more generic like "FIP signal utils"

[alexpiet]

I'll approve it to reduce red tape, but I would clean up the comment I referenced.

[alexpiet]

I would remove mention of bwnm_signal_utils and the benchmark repo

[ksseverson57]

Can do

[ksseverson57]

Do we need to add Johannes as a member of this repo to add him as a reviewer? I don't see him in the drop-down list

[alexpiet]

Do we need to add Johannes as a member of this repo to add him as a reviewer? I don't see him in the drop-down list

Added @j-friedrich to the repo, and requested a review.

@j-friedrich Please review the method developed here by Kyle to estimate signal-to-noise ratio by fitting a moving baseline.

[j-friedrich]

Conceptually, why is an envelope calculated at all? Bleaching (and brightening) is already captured by F0 baseline calculation, and isn't tonic activity signal too? Or is this method supposed to also directly apply to raw traces? The envelope also tends to overfit, see Fig. 4 in your capsule.

Estimating the noise floor robustly from the below-median half of the residual is still biased. Please compare to using the mode instead (using the half sample method), which is less biased. The way the noise is calculated in method deriv is a bit of a strawman to compare to, at the very least use MAD not std. Please add a comparison to the noise estimators we actually already use in aind-ophys-utils.

A default threshold of 1.5 standard deviations seems low — for Gaussian noise, 6.7% of values are above that threshold, though introducing the "rise-rate gate" seems to save you from false positives.

Does the method yield visually reasonable results on real data? Can you run it on a couple of actual FIP traces? If that function is supposed to also work on traces collected with 2-photon imaging, Jinho has assembled a lot of ΔF/F traces in this asset: https://codeocean.allenneuraldynamics.org/data-assets/d9193138-eca9-49a0-b8c8-7975ffcd2faf/test-dataset-for-dff_multiplane-ophys_02

[ksseverson57]

Thanks, Johannes, for the careful review.

Conceptually, the intention of the lower envelope fit is to decompose the dF/F signal into tonic changes on a shorter timescale (~seconds) and a phasic residual, under the assumptions that long-timescale baseline changes are removed and remaining noise is relatively high frequency.

The previous steps in the pipeline should have already removed baseline changes due to bleaching and brightening, as well as motion. My mistake, I should not have mentioned bleaching in top PR message. This method is only intended to be applied to post-baseline corrected dF/F traces downstream of the baseline/motion correction pipeline, not raw F.

Slow tonic changes fit by the envelope can be considered signal. Currently, the SNR measurement only applies to the phasic component. The output could be modified to give separate phasic and tonic SNR estimates and/or added together to estimate total SNR.

I will compare the different noise estimation methods in aind-ophys-utils, test MAD in the simulation. Derivative is a strawman, but it's the current finicky method to which we must compare. Derivative could be saved if it was modified to better estimate noise and more robustly detect peaks (e.g., hyperparameter tuning prominence) so as to not output NaN.

Ultimately, what we need is an interpretable metric that clearly tells us whether FIP signals are above noise level, as well as a graded quality level for day-to-day QC.

I'll attach FIP data assets to the capsule to compare SNR measures on real data. Currently, it's only intended to apply to fiber photometry data, but it could also be applied to averaged ROIs for 2P.

[ksseverson57]

snr_envelope_rr.py: MAD noise estimator; arPLS tonic estimator; total SNR

updated capsule: Kyle_FP_SNR_benchmark_v3 | Code Ocean
https://codeocean.allenneuraldynamics.org/capsule/8413511/tree

Summary

I implemented Johannes' suggestions. Three changes are bundled together
because they touch on separate modules: noise, tonic, and SNR estimation.

1. Default noise estimator changed from folded IQR to MAD. By
   default, use noise_method='aind_mad' from aind-ophys-utils'
   noise_std(method='mad').
2. Default tonic tracker changed from 'envelope' to 'arpls',
   matching what the benchmark notebook found and adopted (envelope
   collapses at high SNR; arPLS doesn't — see notebook Sections 4-5).
3. New: total/tonic/phasic SNR breakdown. result.snr is now
   snr_tonic + snr_phasic rather than phasic-only; tonic and phasic
   components are exposed individually in estimate_components().

Major changes

What	Before	After
noise_method naming	'folded_iqr' (median-anchored)	'aind_mad' (noise_std)
Default tonic_method	'envelope' (valley + spline interp)	'arpls'
result.snr meaning	phasic-only SNR	total SNR (snr_tonic + snr_phasic)

New features

• noise_method='mad_iqr_avg' and mode-anchored 'folded_iqr'.
• tonic_method='arpls' (Asymmetrically Reweighted Penalized Least
  Squares, Baek et al. 2015) — resolves fixed-p ALS's sharp
  low/high-SNR trade-off by re-deriving asymmetry weights from the
  residual's own noise statistics each iteration instead of a single
  fixed p. arPLS outperforms the previous 'envelope' method across
  the tested SNR range.
• Convenience properties estimator.snr_tonic_ / estimator.snr_phasic_
  added.

Not done in this PR (flagging, not fixing)

• No benchmark sweep exists yet for total SNR against ground truth —
  snr_corrected (total) is a derived combination (validated
  snr_phasic_corrected + raw snr_tonic), not an independently-fit
  correction.
• tonic_method='envelope'/'als' remain implemented and directly
  comparable but are no longer the default; no changes to their
  behavior in this PR.
• Application on real data. Need a test on FIP data using various
  sensors and regions, should show consistent values within subject.
  Could also do real data + tonic and/or phasic noise injection and
  test scaling.