Enable IRASA multi-method support and add regression tests

[camiloberutti]

Summary

This PR fixes compute_irasa in neurodsp.aperiodic.irasa to support all spectral estimation methods (welch, medfilt, multitaper, wavelet) instead of only welch. Currently, compute_irasa unconditionally injects the nperseg parameter, which is invalid for non-Welch methods, causing AssertionError or KeyError.

Problem Statement

Users cannot use compute_irasa with methods other than Welch because:

1. Missing method-aware parameter injection: The function always adds nperseg regardless of the spectral method, violating parameter constraints for medfilt and multitaper.
2. Missing multitaper entry in SPECTRUM_INPUTS: The parameter validation dictionary in power.py lacks a multitaper key, causing KeyError when the method is used.
3. Frequency-grid misalignment: Different spectral methods produce output arrays of different sizes. IRASA resampling needs robust frequency-grid alignment.

Error Examples (Before Patch)

# Medfilt: AssertionError
>>> compute_irasa(sig, fs, method='medfilt', filt_len=0.0001)
AssertionError: Parameter nperseg not expected for medfilt estimation method

# Multitaper: KeyError
>>> compute_irasa(sig, fs, method='multitaper', bandwidth=1.0)
KeyError: 'multitaper'

Changes

File 1: neurodsp/aperiodic/irasa.py

Change 1.1: Method-aware nperseg injection (line ~72)

# BEFORE
if 'nperseg' not in spectrum_kwargs:
    spectrum_kwargs['nperseg'] = int(4 * fs)

# AFTER
# Only Welch uses `nperseg`; avoid injecting it for non-Welch methods.
if spectrum_kwargs.get('method', 'welch') == 'welch' and 'nperseg' not in spectrum_kwargs:
    spectrum_kwargs['nperseg'] = int(4 * fs)

Change 1.2: Robust frequency-grid alignment for resampled spectra (lines ~86–96)

# BEFORE
freqs_up, psd_up = compute_spectrum(sig_up, h_val * fs, **spectrum_kwargs)
freqs_dn, psd_dn = compute_spectrum(sig_dn, fs / h_val, **spectrum_kwargs)
psds[ind, :] = np.sqrt(psd_up * psd_dn)  # assumes same length

# AFTER
freqs_up, psd_up = compute_spectrum(sig_up, h_val * fs, **spectrum_kwargs)
freqs_dn, psd_dn = compute_spectrum(sig_dn, fs / h_val, **spectrum_kwargs)

# Align spectra to the original frequency grid for methods whose output length
# changes with signal length (for example medfilt and multitaper).
psd_up_i = np.interp(freqs, freqs_up, psd_up, left=np.nan, right=np.nan)
psd_dn_i = np.interp(freqs, freqs_dn, psd_dn, left=np.nan, right=np.nan)

# Calculate the geometric mean of h and 1/h on a shared frequency grid.
psds[ind, :] = np.sqrt(psd_up_i * psd_dn_i)

Change 1.3: Use np.nanmedian instead of np.median (line ~104)

# BEFORE
psd_aperiodic = np.median(psds, axis=0)

# AFTER
psd_aperiodic = np.nanmedian(psds, axis=0)

File 2: neurodsp/spectral/power.py

Change 2.1: Add multitaper to SPECTRUM_INPUTS dictionary (line ~80)

# BEFORE
SPECTRUM_INPUTS = {
    'wavelet' : ['freqs', 'avg_type', 'n_cycles', 'scaling', 'norm'],
    'fft' : ['window', 'f_range'],
    'welch' : ['avg_type', 'window', 'nperseg', 'noverlap', 'nfft', 'fast_len', 'f_range'],
    'medfilt' : ['filt_len', 'f_range'],
}

# AFTER
SPECTRUM_INPUTS = {
    'wavelet' : ['freqs', 'avg_type', 'n_cycles', 'scaling', 'norm'],
    'fft' : ['window', 'f_range'],
    'welch' : ['avg_type', 'window', 'nperseg', 'noverlap', 'nfft', 'fast_len', 'f_range'],
    'medfilt' : ['filt_len', 'f_range'],
    'multitaper' : ['bandwidth', 'n_tapers', 'low_bias', 'eigenvalue_weighting'],
}

Testing

New test file: test_irasa_multimethod.py (see attached)

• Tests compute_irasa with all four methods
• Verifies output consistency
• Validates backward compatibility

Usage example after patch:

from neurodsp.aperiodic import compute_irasa
import numpy as np

sig = np.random.randn(86400)
fs = 1.0
f_range = [1e-5, 1e-3]

# All now work
freqs_w, ap_w, pe_w = compute_irasa(sig, fs, f_range=f_range, method='welch', nperseg=2048)
freqs_m, ap_m, pe_m = compute_irasa(sig, fs, f_range=f_range, method='medfilt', filt_len=0.0001)
freqs_t, ap_t, pe_t = compute_irasa(sig, fs, f_range=f_range, method='multitaper', bandwidth=2.0)
freqs_w, ap_w, pe_w = compute_irasa(sig, fs, f_range=f_range, method='wavelet', freqs=np.logspace(-5, -3, 40))