Hi MNE team,
I am currently trying to create an automated pipeline that requires the cardiac field artefact to be removed from the R-peak-based epochs. I have been using find_bads_ecg on 1000 Hz EEG and noticed it flags different components depending on whether I resample first. I was planning on using mostly the defaults so method='ctps' and threshold='auto' are what I was going to use.

_get_ctps_threshold seems to be using the sampling rate as N (N = self.info["sfreq"]) in Eq. (13) from the CTPS paper that the documentation refers to (https://doi.org/10.1109/TBME.2008.926677).

This is indeed in the documentation itself ("will compute the threshold automatically based on the sampling frequency"), so it looks deliberate. But I don't think N is meant to be the sampling rate. The paper's own wording, "N being the number of data points", is ambiguous read on its own, but the explanation before it isn't:

"Cross trial distributions are time-dependent histograms (in our case, phase histograms) that are calculated across trials (data segments) for each time point relative to the onset of an event (e.g., the R-peak)."

At each time point the distribution under test has one sample per trial, so N in Eq. (13) should be (if I understand it correctly) the number of trials. Time is the axis the test is repeated over, not the one it's computed across. The kuiper() function in ctps_.py reads it the same way, using n_trials = shape[0].

The effect is that the threshold is based on the how the data was recorded and not on the number of heartbeat epochs that are used in the calculation:

• 250 Hz -> 0.32 threshold
• 500 Hz -> 0.23 threshold
• 1000 Hz -> 0.16 threshold

The scores don't move with sampling rate, so a component that is clearly locked to the R-peak can pass at 1000 Hz and fail at 250 Hz on the same data. I can see that sampling rate can have an effect though, since the scores are a max over time points and the null creeps up a little at higher rates, but that would look more like a higher threshold at higher rates rather than a lower one.

Once seeing this, I used claude (AI disclaimer for the following part) to see if anything else looked suspicious with the ctps implementation and it mentioned that there could also be a discrepancy with the scores themselves. More specifically that the scores from _prob_kuiper are normalized to [0, 1] and tend to V**2, whereas _get_ctps_threshold returns a raw Kuiper statistic V, so the comparison in find_bads_ecg also looks like "it's off by a square". Changing only N leaves that in place and ends up stricter than the current behaviour, not closer to the paper. Lastly it mentioned that "kuiper() already returns the raw statistic alongside the normalized one (ks_dynamics), but find_bads_ecg discards it and keeps only pk_dynamics."

Finally, do note that I am unsure if I am just misreading the implementation or the original paper (and if the AI is too). But please let me know if I can provide anything else that could be useful.