Fix psd_array_welch for good data spans shorter than n_per_seg (#13039)

[CedricConday]

Reference issue

Fixes #13039.

What does this implement/fix?

When bad_* annotations split the data into good spans, a span shorter than n_per_seg previously raised from SciPy:

ValueError: noverlap must be less than nperseg.

Such spans are now analysed with a window shrunk to the span's own length, with n_overlap clamped below it and a warning about the reduced resolution. n_fft is left unchanged, so every span is zero-padded to the same FFT length and shares one frequency grid — short spans simply get coarser spectral resolution. No data is discarded.

An explicit ndarray window cannot be shortened to match, so that combination raises a clear ValueError pointing at the remedy (pass a shorter window, or reduce n_per_seg/n_fft) rather than letting SciPy raise a length mismatch.

Why shrinking rather than zero-padding

An earlier revision of this PR zero-padded short spans up to n_per_seg instead, to match the direction of scipy/scipy#25608. @viranovskaya's simulations showed that padding is not neutral: it underestimates alpha power by 5.8% / 13.7% / 21.5% as the share of short spans rises to 20% / 50% / 80%, and the SciPy PR's right-padding reproduces that bias exactly (matched to 2.2e-16). Shrinking stayed within 0.2%.

Padding a good span also asserts something untrue about the data — that the signal was zero there — when in fact those samples were excluded as bad. Shrinking asserts only that fewer samples are available, so the resolution is coarser there, which is both true and visible to the user via the warning.

Test

• test_psd_welch_short_span_kept — the warning-and-keep path (one short span among long ones) and the all-spans-short case.
• test_psd_welch_short_span_array_window_raises — the fixed-length window array case.
• test_psd_welch_short_span_recovers_band_power — asserts the recovered alpha-band power matches the uninterrupted signal within 10%. This is the regression guard for the bias above: it passes at −0.3% with shrinking and fails at −78% against the zero-padding implementation, which a finiteness check cannot distinguish.

Full mne/time_frequency/tests/test_psd.py green (23 passed).

---

Disclosure: I'm an AI engineer and use Claude Code in my workflow. All code here is hand-reviewed and the reasoning is mine; I'm happy to walk through any of it.

[welcome]

Hello! 👋 Thanks for opening your first pull request here! ❤️ We will try to get back to you soon. 🚴

[CarinaFo]

Hi Cedric, thank you for your contribution.

Could you disclose your AI usage according to our adopted policy.

I ran into this same issue in my analysis recently and ended up deciding that dropping short spans was the safer approach, so I would be rather cautious about shrinking noverlap, particularly without warning the user about this behaviour.
I suggest we sync with @drammock on the right approach here (I haven't looked at the "old" behaviour but I think it would be important to know why noverlap < nperseg did not throw a ValueError in older MNE versions.)

[CedricConday]

Hi Carina, thanks for the careful review.

On the AI disclosure — I've added it to the PR description per the policy: I'm an AI engineer and use Claude Code in my workflow; I find, fix, and test, then review and verify everything before opening it under my name. Happy to walk through any part of the reasoning.

On the approach — I agree, silently shrinking noverlap didn't sit right with me either. I've pushed an update that instead drops good-data spans shorter than n_per_seg with a warning, and raises a clear ValueError if every span is too short. Glad that matches where you landed.

Happy to sync with @drammock before this is finalized. On your open question — I don't want to guess at why noverlap < n_per_seg didn't raise in older MNE, so I'll dig into the history and report back here, so we're deciding on the real prior behavior rather than assumptions.

Thanks!

[CedricConday]

Quick follow-up on the history question — I traced it. In 2015 (ca23ab6) _check_nfft silently clamped both: oversized n_fft down to the data length, and n_overlap >= n_fft down to n_fft - 1 — nothing raised. The ValueError for n_overlap >= n_per_seg was added in 2017 with #4003 (the PR that introduced n_per_seg). That same PR kept the segment-length clamp silent — n_per_seg = n if n_per_seg > n — and that line is still on main. So the short-span case never raised: the window just gets quietly shrunk to the span length, which is exactly the behavior here. Only the overlap check ever became strict; the segment-length clamp stayed silent throughout.

[drammock]

So if I'm understanding correctly, the problem is:

• we effectively hard-code the values of n_overlap, n_per_seg, and n_fft that we pass to scipy.signal.spectrogram, by baking them into a functools.partial, here:

  mne-python/mne/time_frequency/psd.py

  Lines 247 to 256 in 31a7326

  	_func = partial(
  	spectrogram,
  	detrend=detrend,
  	noverlap=n_overlap,
  	nperseg=n_per_seg,
  	nfft=n_fft,
  	fs=sfreq,
  	window=window,
  	mode=mode,
  	)

• we do that before we've sorted out which spans of the signal need to be skipped / which spans we're going to analyze (AKA x_splits, here and here)
• so if one of the spans is shorter than n_per_seg, scipy will change n_per_seg internally but will not change n_overlap, and then it will error out a bit further on

I lean toward considering this a SciPy bug: the user passed input that ought to work, and then scipy-internal code changed one value, then chokes because that value is no longer OK relative the other passed value. Would love @larsoner's opinion as to whether upstreaming this report makes sense.

As for what to do about it on our end (mostly repeating what others have said here), we could:

1. error out (current behavior)
2. drop short segments (with a warning)
3. catch this problem before passing a span to scipy.signal.spectrogram and adjust the values of both n_per_seg and n_overlap in the partial (workaround)
4. report upstream and (wait for a) fix; this is a severe enough bug though that even if this is the long term approach I think we should also do (3) in the meantime.

I lean toward (3) or (4 & 3)

[CedricConday]

Thanks @drammock — your root-cause read is exactly right: n_per_seg/n_overlap/n_fft get baked into the functools.partial before spans are split, so a span shorter than n_per_seg makes SciPy clamp nperseg internally but leave noverlap, then error.

I've gone with option 3. Instead of dropping short spans, each span now gets a per-span partial that overrides the baked-in values: nperseg shrunk to the span length and noverlap clamped below it. n_fft is untouched, so every span is zero-padded to the same length and shares one frequency grid — short spans just get coarser spectral resolution, which I surface with a warning. No data is discarded. The all-long path is unchanged (the test_bad_annot_handling equivalence test still matches at rtol=1e-15).

Two things worth your call:

1. Array windows. For string/tuple windows SciPy regenerates the window at the shrunk nperseg, so this is transparent. An explicit ndarray window has a fixed length and can't be shrunk — that combination (array window + a span shorter than n_per_seg) isn't handled here. Prefer I drop-with-warning in just that sub-case, or raise a clearer error?
2. Upstream (option 4). Happy to open a SciPy issue for the clamp-nperseg-but-not-noverlap behavior. Want me to, and link it here?

[drammock]

Happy to open a SciPy issue for the clamp-nperseg-but-not-noverlap behavior. Want me to, and link it here?

Yeah I think it's worth asking them if they consider it bug-ish and want it fixed.

An explicit ndarray window has a fixed length and can't be shrunk

I don't have a clear sense of what users would prefer here (it might vary). drop-segment-plus-warning is a bit more friendly I guess, but would still be catastrophic if a particular dataset was mostly/all short spans... but I guess in either case the remedy is "pass in a shorter custom window array" so maybe that's fine.

[CedricConday]

Pushed ae6fcbf handling the array-window case.

Array/ndarray window + short span — went with the clear-error route. A named/tuple window is regenerated by SciPy at the shrunk nperseg, so those spans stay transparent; but a fixed-length window array can't be shortened, so instead of letting SciPy raise a cryptic length-mismatch, psd_array_welch now raises an actionable ValueError pointing the user to pass a shorter window array or reduce n_per_seg/n_fft. I preferred this over drop-with-warning since silently dropping is catastrophic on a mostly-short dataset (your concern), and the remedy you named — "pass a shorter custom window array" — is exactly what the error tells them to do. Regression test added; full test_psd.py green locally (22/22).

Upstream (option 4) — I'll open a SciPy issue for the clamp-nperseg-but-not-noverlap behavior and link it here.

[CedricConday]

Opened the upstream SciPy issue: scipy/scipy#25608 — clamp-nperseg-but-not-noverlap for short input, with a minimal reproducer. Our per-span workaround here stands regardless of how they decide to handle it upstream.

[CarinaFo]

Nice work, the implementation looks good @CedricConday. I had a look at the scipy issue and it seems the suggested fix from scipy uses a different approach: the original signal is zero-padded up to nperseg (with nperseg left unchanged), rather than shrinking nperseg/window to the signal length as this PR does (The scipy PR isn't reviewed yet, so this could still change.)
I guess we have to decide on a strategy (zero padding vs. shrinking) and consider whether this might be problem for the user if our preferred solution does not align with scipy's. Worth further discussion.

[viranovskaya]

To make this choice less abstract, I ran a small simulation comparing the three options.

I generated 200 stationary signals with known 10 Hz and 20 Hz activity, split them into long and short good-data spans, and compared each result with the PSD of the uninterrupted signal.

In this test:

• shrinking the window was closest to the original PSD and kept the median alpha and beta errors below 0.2%;
• zero-padding increasingly underestimated alpha and beta power: by about 6% when 20% of the samples were in short spans, 14% at 50%, and 21% at 80%;
• dropping gave reasonable bandpower estimates in this stationary example, but discarded more data as short spans became common.

This does not mean that shrinking is always the best choice. Real short spans may contain a different signal, and that would need a separate test. But the result suggests that zero-padding is not neutral.

It may be useful to add a regression test that checks the estimated power of a known signal, rather than only checking that the output is finite. The PR description also still says that short spans are dropped, while the current code zero-pads them.

I can share the script and plot if useful.

[CarinaFo]

Hi @viranovskaya,

would be really useful for us if you could share the code and plots for the simulations so we can make an informed decision on what is the better approach here.

Thank your for your great work.

Carina

[viranovskaya]

Hi Carina,

Thank you. I’ve attached the plot and a ZIP containing the simulation script, full CSV results, and a short summary.

In this stationary simulation, shrinking produced the lowest full-spectrum error. Zero-padding increasingly underestimated band power as short spans became more common, while dropping preserved band power reasonably well but discarded more data.

This test does not show that shrinking is always preferable, since real short spans may contain a different signal. It does show that zero-padding is not neutral.

[drammock]

I’ve attached the plot and a ZIP containing the simulation script, full CSV results, and a short summary

zip file attachment didn't succeed; can you post the script and CSV as as a GitHub Gist maybe? Or just the script is enough, if it regenerates the CSV.

[viranovskaya]

Thanks — the ZIP upload failed. I’ve posted the full script as a public Gist here:

https://gist.github.com/viranovskaya/f245abb502567c2dbc665b59daad7284

Running python run_simulation.py regenerates simulation_results.csv, summary.md, and strategy_comparison.png in the same directory. It requires NumPy, SciPy, and Matplotlib.

[viranovskaya]

Thanks — I updated the simulation to separate end-padding, centred padding, and centred padding with a valid-window correction.

The result changes the earlier interpretation. With 80% of samples in short spans, median alpha error was -21.5% for end-padding, -6.5% for centred padding, and -0.3% after the correction. The corrected centred condition also had lower median log-PSD RMSE than shrinking in this stationary test (0.0434 vs 0.0585), although median beta error was +1.3%.

So the bias is not a property of zero-padding alone; placement and normalisation both matter. This is still a stationary synthetic test and does not decide what is best when short spans contain a different physiological state.

The updated script and full table are in the same Gist:
https://gist.github.com/viranovskaya/f245abb502567c2dbc665b59daad7284

[CarinaFo]

Thank you for the detailed simulations, this does change the picture quite a lot. I would still opt for shrinking because it is less complex but then scipy does zero pad, so would be also nice to be consistent. Next up would be to compare our way of zero padding (centered and corrected) with the scipy PR

[viranovskaya]

Thanks — I ran the current SciPy PR implementation directly on the same 200 simulations.

The PR right-pads each short input to nperseg, and its output matched my earlier end-padding condition (maximum absolute PSD difference: 2.2e-16).

At 20%, 50%, and 80% short samples, median alpha error was −5.8%, −13.7%, and −21.5%, respectively. Shrinking remained close to zero, while centred padding with the valid-window correction was also close to zero and had lower log-PSD RMSE in this stationary simulation.

So the SciPy PR would give consistent behaviour, but the right-padding itself is not neutral here. I’ll add the direct comparison script, table, and plot to the Gist.

[CedricConday]

Sorry for the silence — picking this back up.

First, a correction to my own PR that I should have flagged sooner: the branch currently end-pads, and @viranovskaya's numbers say that is the worst of the three options. 4ebb31b65 pads with np.pad(span, [(0, n_per_seg - span_len)]), which is exactly the end-padding condition she matched to the SciPy PR at 2.2e-16. So the code as it stands carries the −5.8% / −13.7% / −21.5% alpha bias. That needs to change regardless of which direction we pick.

My recommendation: go back to shrinking.

@CarinaFo, I think the consistency argument dissolves once we know what SciPy's approach actually costs. Consistency is only worth having when the thing we are consistent with is unbiased, and viranovskaya's measurement says right-padding is not: it underestimates alpha power by up to 21.5% in a stationary test. Matching SciPy here would mean adopting a known bias in order to agree with it.

There is also a reason the right answer should differ from SciPy's, beyond the numbers. SciPy is handed a single short array and has to decide what that array means; padding it is a defensible reading. We are in a different situation: our short spans are carved out of a longer recording by bad_* annotations. Zero-padding one asserts that the signal was zero over the padded region — but it was not zero, it was excluded as bad. That is a false statement about the data, and it injects exactly the power deficit that was measured. Shrinking asserts something true instead: we have fewer samples here, so we get coarser spectral resolution here. Same reason we would not pad across a bad segment in the middle of a span.

The failure modes also differ in kind. Shrinking's cost is coarser resolution on short spans — visible, warned about, and interpretable. Padding's cost is a silent underestimate of band power, which is invisible and lands on precisely the quantity most people compute from this function. Given the choice, I would rather be transparently coarse than quietly wrong.

And it happens to be the simplest option, which was your instinct anyway — the usual accuracy/complexity trade-off doesn't bite here. n_fft stays fixed, so all spans still share one frequency grid; only the window length varies.

The one option I would not dismiss is centred padding with the valid-window correction, which viranovskaya measured at −0.3% alpha and a slightly better log-PSD RMSE than shrinking. If we want maximum accuracy I would take that over end-padding without hesitation. I lean against it only because it is a third behaviour — consistent with neither SciPy nor the simple approach — and it adds both a placement rule and a normalisation correction to maintain, for a band-power gain that is within noise of shrinking in this test.

Two follow-ups either way:

1. @viranovskaya's suggestion to assert recovered power of a known signal rather than just finiteness is the right call, and it is the test that would have caught this. I'll add a band-power regression test with the chosen strategy.
2. The PR description is still describing the drop behaviour from two implementations ago. I'll rewrite it to match whatever we land on.

@viranovskaya — thank you, genuinely. Running the SciPy PR against the same 200 simulations is what turned this from a preference argument into a decidable one, and it caught a regression in my own branch.

Happy to push the change as soon as you and @drammock agree on the direction.

AI-assisted, human-reviewed.

[CedricConday]

Pushed the decision: back to shrinking, plus the band-power regression test @viranovskaya asked for.

3d3f15422 reverts the zero-padding and restores the per-span shrunk window (and with it the clear ValueError for a fixed-length window array, which padding had made unnecessary).

The new test_psd_welch_short_span_recovers_band_power is the guard that was missing. It chops a known 10 Hz signal into good spans all shorter than n_per_seg and asserts the recovered alpha-band power matches the uninterrupted estimate within 10%. Run both ways on identical input:

• shrinking: −0.26% → passes
• zero-padding: −77.8% → fails (Max relative difference: 0.778)

That gap is the point — both produce finite, plausible-looking spectra, so the old finiteness assertions could not tell them apart.

I've also rewritten the PR description, which was still describing the drop behaviour from two implementations ago.

@CarinaFo — this lands where you leaned originally. I don't think the SciPy-consistency argument survives the measurement: matching them here means importing a known bias, and our situation differs from theirs anyway, since our short spans are carved out of a longer recording rather than being the whole input.

Full test_psd.py green (23 passed).

AI-assisted, human-reviewed.