From 0af5d36b8d359834e2f23c2db97269ada6531c0c Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Sun, 12 Oct 2025 15:01:09 +0100
Subject: [PATCH 1/8] add compute_spectrum_fft func

---
 neurodsp/spectral/power.py            | 39 +++++++++++++++++++++++++--
 neurodsp/spectral/utils.py            | 32 ++++++++++++++++++++++
 neurodsp/tests/spectral/test_power.py |  5 ++++
 neurodsp/tests/spectral/test_utils.py | 21 +++++++++++++++
 4 files changed, 95 insertions(+), 2 deletions(-)

diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py
index 07284ff45..ffc23d6de 100644
--- a/neurodsp/spectral/power.py
+++ b/neurodsp/spectral/power.py
@@ -16,7 +16,7 @@
 from neurodsp.utils.checks import check_param_options
 from neurodsp.utils.outliers import discard_outliers
 from neurodsp.timefrequency.wavelets import compute_wavelet_transform
-from neurodsp.spectral.utils import trim_spectrum, window_pad
+from neurodsp.spectral.utils import trim_spectrum, window_pad, get_positive_fft_outputs
 from neurodsp.spectral.checks import check_spg_settings, check_mt_settings
 
 ###################################################################################################
@@ -71,9 +71,10 @@ def compute_spectrum(sig, fs, method='welch', **kwargs):
 
 
 SPECTRUM_INPUTS = {
+    'wavelet' : ['freqs', 'avg_type', 'n_cycles', 'scaling', 'norm'],
+    'fft' : ['f_range'],
     'welch' : ['avg_type', 'window', 'nperseg', 'noverlap', 'nfft', \
                'fast_len', 'f_range', 'outlier_percent'],
-    'wavelet' : ['freqs', 'avg_type', 'n_cycles', 'scaling', 'norm'],
     'medfilt' : ['filt_len', 'f_range'],
 }
 
@@ -137,6 +138,40 @@ def compute_spectrum_wavelet(sig, fs, freqs, avg_type='mean', **kwargs):
     return freqs, spectrum
 
 
+@multidim(select=[0])
+def compute_spectrum_fft(sig, fs, f_range=None):
+    """Compute the power spectrum based on a single FFT.
+
+    Parameters
+    ----------
+    sig : array
+        Time series.
+    fs : float
+        Sampling rate, in Hz.
+    f_range : list of [float, float], optional
+        Frequency range to sub-select from the power spectrum.
+
+    Returns
+    -------
+    freqs : 1d array
+        Frequencies at which the measure was calculated.
+    spectrum : array
+        Power spectral density.
+    """
+
+    # Compute the FFT and take the real part of the FFT power estimate
+    spectrum = np.real(np.fft.fft(sig))
+
+    # Compute the frequency vector, and extract positive frequency & power values
+    freqs = np.fft.fftfreq(len(sig), 1/fs)
+    freqs, spectrum = get_positive_fft_outputs(freqs, spectrum)
+
+    if f_range:
+        freqs, spectrum = trim_spectrum(freqs, spectrum, f_range)
+
+    return freqs, spectrum
+
+
 def compute_spectrum_welch(sig, fs, avg_type='mean', window='hann',
                            nperseg=None, noverlap=None, nfft=None,
                            fast_len=False, f_range=None, outlier_percent=None):
diff --git a/neurodsp/spectral/utils.py b/neurodsp/spectral/utils.py
index b32614bc1..1aa62d8f3 100644
--- a/neurodsp/spectral/utils.py
+++ b/neurodsp/spectral/utils.py
@@ -130,6 +130,38 @@ def trim_spectrogram(freqs, times, spg, f_range=None, t_range=None):
     return freqs_ext, times_ext, spg_ext
 
 
+def get_positive_fft_outputs(freqs, powers=None, drop_zero=False):
+    """Get the positive frequency values for an FFT.
+
+    Parameters
+    ----------
+    freqs : 1d array
+        Frequency vector corresponding to the FFT estimate, with positive & negative frequencies.
+    powers : 1d array, optional
+        Complex power value estimates from the FFT.
+    drop_zero : bool, optional, default: False
+        Whether to drop the estimate for frequency of 0.
+
+    Returns
+    -------
+    freqs : 1d array
+        Frequencies at which the measure was calculated.
+    spectrum : array
+        Power spectral density.
+        Only returned if an input power spectrum is passed.
+    """
+
+    start_ind = 1 if drop_zero else 0
+
+    # Get the max positive ind as half length, rounded up if the length is odd
+    end_ind = int(np.ceil(len(freqs) / 2))
+
+    if powers is not None:
+        return freqs[start_ind:end_ind], powers[start_ind:end_ind]
+    else:
+        return freqs[start_ind:end_ind]
+
+
 def window_pad(sig, nperseg, noverlap, npad, fast_len,
                nwindows=None, nsamples=None, pad_left=None, pad_right=None):
     """Pads windows (for Welch's PSD) with zeros.
diff --git a/neurodsp/tests/spectral/test_power.py b/neurodsp/tests/spectral/test_power.py
index 637d0ed82..849e40182 100644
--- a/neurodsp/tests/spectral/test_power.py
+++ b/neurodsp/tests/spectral/test_power.py
@@ -55,6 +55,11 @@ def test_compute_spectrum_2d(tsig2d):
     assert freqs.shape[-1] == spectrum.shape[-1]
     assert spectrum.ndim == 2
 
+def test_compute_spectrum_fft(tsig, tsig_sine):
+
+    freqs, spectrum = compute_spectrum_fft(tsig, FS)
+    assert freqs.shape == spectrum.shape
+
 def test_compute_spectrum_welch(tsig, tsig_sine):
 
     freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='mean')
diff --git a/neurodsp/tests/spectral/test_utils.py b/neurodsp/tests/spectral/test_utils.py
index 430c90e09..7a6ec0e10 100644
--- a/neurodsp/tests/spectral/test_utils.py
+++ b/neurodsp/tests/spectral/test_utils.py
@@ -43,6 +43,27 @@ def test_trim_spectrogram():
     assert_equal(f_ext, np.array([6, 7, 8]))
     assert_equal(t_ext, times)
 
+def test_get_positive_fft_outputs():
+
+    # Test odd length
+    freqs_odd = np.array([0, 1, 2, -2, -1])
+    powers_odd = np.array([0, 1, 2, 3, 4])
+    freqs_out_odd, powers_out_odd = get_positive_fft_outputs(freqs_odd, powers_odd)
+    assert freqs_out_odd.shape == powers_out_odd.shape
+    assert np.array_equal(freqs_out_odd, freqs_odd[:3])
+
+    # Test drop zero
+    freqs_out_odd2, powers_out_odd2 = get_positive_fft_outputs(freqs_odd, powers_odd, drop_zero=True)
+    assert freqs_out_odd2[0] != 0
+    assert powers_out_odd2[0] != 0
+
+    # Test even length
+    freqs_even = np.array([0, 1, 2, -3, -2, -1])
+    powers_even = np.array([0, 1, 2, 3, 4, 5])
+    freqs_out_even, powers_out_even = get_positive_fft_outputs(freqs_odd, powers_odd)
+    assert freqs_out_even.shape == powers_out_even.shape
+    assert freqs_out_even.shape == powers_out_even.shape
+    assert np.array_equal(freqs_out_even, freqs_even[:3])
 
 @pytest.mark.parametrize("fast_len", [True, False])
 def test_window_pad(fast_len):

From 3ae0c6879bb330201490285f677a8748db52fb6f Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Sun, 12 Oct 2025 15:09:22 +0100
Subject: [PATCH 2/8] use compute_spectrum_fft in medfilt

---
 neurodsp/spectral/power.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py
index ffc23d6de..b716bd3af 100644
--- a/neurodsp/spectral/power.py
+++ b/neurodsp/spectral/power.py
@@ -159,11 +159,11 @@ def compute_spectrum_fft(sig, fs, f_range=None):
         Power spectral density.
     """
 
-    # Compute the FFT and take the real part of the FFT power estimate
-    spectrum = np.real(np.fft.fft(sig))
+    # Compute the FFT and compute power
+    spectrum = np.abs(np.fft.fft(sig))**2.
 
     # Compute the frequency vector, and extract positive frequency & power values
-    freqs = np.fft.fftfreq(len(sig), 1/fs)
+    freqs = np.fft.fftfreq(len(sig), 1. / fs)
     freqs, spectrum = get_positive_fft_outputs(freqs, spectrum)
 
     if f_range:
@@ -299,16 +299,15 @@ def compute_spectrum_medfilt(sig, fs, filt_len=1., f_range=None):
     >>> freqs, spec = compute_spectrum_medfilt(sig, fs=500)
     """
 
-    # Take the positive half of the spectrum, since it's symmetrical
-    ft = np.fft.fft(sig)[:int(np.ceil(len(sig) / 2.))]
-    freqs = np.fft.fftfreq(len(sig), 1. / fs)[:int(np.ceil(len(sig) / 2.))]
+    # Compute spectrum estimate as a single FFT
+    freqs, spectrum = compute_spectrum_fft(sig, fs)
 
     # Convert median filter length from Hz to samples, and make sure it is odd
     filt_len_samp = int(filt_len / (freqs[1] - freqs[0]))
     if filt_len_samp % 2 == 0:
         filt_len_samp += 1
 
-    spectrum = medfilt(np.abs(ft)**2. / (fs * len(sig)), filt_len_samp)
+    spectrum = medfilt(spectrum / (fs * len(sig)), filt_len_samp)
 
     if f_range:
         freqs, spectrum = trim_spectrum(freqs, spectrum, f_range)

From 7a0b6b255ad5538239d20442bad28c62b672d350 Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Sun, 12 Oct 2025 15:56:47 +0100
Subject: [PATCH 3/8] use rfft for compute_spectrum_fft & lints

---
 neurodsp/spectral/power.py            | 11 ++++-------
 neurodsp/spectral/utils.py            | 19 ++++++++++---------
 neurodsp/tests/spectral/test_power.py | 15 +++++----------
 3 files changed, 19 insertions(+), 26 deletions(-)

diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py
index b716bd3af..972992a64 100644
--- a/neurodsp/spectral/power.py
+++ b/neurodsp/spectral/power.py
@@ -159,12 +159,9 @@ def compute_spectrum_fft(sig, fs, f_range=None):
         Power spectral density.
     """
 
-    # Compute the FFT and compute power
-    spectrum = np.abs(np.fft.fft(sig))**2.
-
-    # Compute the frequency vector, and extract positive frequency & power values
-    freqs = np.fft.fftfreq(len(sig), 1. / fs)
-    freqs, spectrum = get_positive_fft_outputs(freqs, spectrum)
+    # Compute the FFT and convert to power & compute corresponding frequency vector
+    spectrum = np.abs(np.fft.rfft(sig)) ** 2.
+    freqs = np.fft.rfftfreq(len(sig), 1. / fs)
 
     if f_range:
         freqs, spectrum = trim_spectrum(freqs, spectrum, f_range)
@@ -374,7 +371,7 @@ def compute_spectrum_multitaper(sig, fs, bandwidth=None, n_tapers=None,
                              "Could not compute spectrum with low_bias=True.")
 
     # Compute Fourier transform on signal weighted by each slepian sequence
-    freqs = np.fft.rfftfreq(sig_len, 1. /fs)
+    freqs = np.fft.rfftfreq(sig_len, 1. / fs)
     spectra = np.abs(np.fft.rfft(slepian_sequences[:, np.newaxis] * sig)) ** 2
 
     # combine estimates to compute final spectrum
diff --git a/neurodsp/spectral/utils.py b/neurodsp/spectral/utils.py
index 1aa62d8f3..051a8e395 100644
--- a/neurodsp/spectral/utils.py
+++ b/neurodsp/spectral/utils.py
@@ -149,6 +149,11 @@ def get_positive_fft_outputs(freqs, powers=None, drop_zero=False):
     spectrum : array
         Power spectral density.
         Only returned if an input power spectrum is passed.
+
+    Notes
+    -----
+    This can be used to extract positive only frequency from an FFT, for example,
+    as returned by `np.fft.fft` & np.fft.fftfreq`.
     """
 
     start_ind = 1 if drop_zero else 0
@@ -176,7 +181,7 @@ def window_pad(sig, nperseg, noverlap, npad, fast_len,
         Number of points to overlap between segments, applied prior to zero padding.
     npad : int
         Number of samples to zero pad windows per side.
-    fast_len : bool, optional
+    fast_len : bool
         Moves nperseg to the fastest length to reduce computation.
         Adjusts zero-padding to account for the new nperseg.
         See scipy.fft.next_fast_len for details.
@@ -190,13 +195,11 @@ def window_pad(sig, nperseg, noverlap, npad, fast_len,
     """
 
     if sig.ndim == 2:
-        # Determine the number of samples and padding once,
-        #   to prevent redundant computation in the loop
+
+        # Determine nsamples & padding once, to prevent redundant computation in the loop
         nwindows = int(np.ceil(len(sig[0])/nperseg))
         if nsamples is None or pad_left is None or pad_right is None:
-            nsamples, pad_left, pad_right = _find_pad_size(
-                nperseg, npad, fast_len
-            )
+            nsamples, pad_left, pad_right = _find_pad_size(nperseg, npad, fast_len)
 
         # Recursively call window_pad on each signal
         for sind, csig in enumerate(sig):
@@ -226,9 +229,7 @@ def window_pad(sig, nperseg, noverlap, npad, fast_len,
 
         if nsamples is None or pad_left is None or pad_right is None:
             # Skipped if called from the 2d case
-            nsamples, pad_left, pad_right = _find_pad_size(
-                nperseg, npad, fast_len
-            )
+            nsamples, pad_left, pad_right = _find_pad_size(nperseg, npad, fast_len)
 
         # Window signal
         sig_windowed = np.zeros((nwindows, nsamples))
diff --git a/neurodsp/tests/spectral/test_power.py b/neurodsp/tests/spectral/test_power.py
index 849e40182..5df0a3342 100644
--- a/neurodsp/tests/spectral/test_power.py
+++ b/neurodsp/tests/spectral/test_power.py
@@ -71,8 +71,7 @@ def test_compute_spectrum_welch(tsig, tsig_sine):
     # Use a rectangular window with a width of one period/cycle and no overlap
     #   The spectrum should just be a dirac spike at the first frequency
     window = np.ones(FS)
-    _, psd_welch = compute_spectrum(tsig_sine, FS, method='welch',
-                                    nperseg=FS, noverlap=0, window=window)
+    _, psd_welch = compute_spectrum_welch(tsig_sine, FS, nperseg=FS, noverlap=0, window=window)
 
     # Spike at frequency 1
     assert np.abs(psd_welch[FREQ_SINE] - 0.5) < EPS
@@ -86,7 +85,7 @@ def test_compute_spectrum_welch(tsig, tsig_sine):
     assert np.allclose(psd_welch[0:FREQ_SINE], expected_answer, atol=EPS)
 
     # Test zero padding
-    freqs, spectrum = compute_spectrum(
+    freqs, spectrum = compute_spectrum_welch(
         np.tile(tsig, (2, 1)), FS, nperseg=100, noverlap=0, nfft=1000, f_range=(1, 200)
     )
     assert np.all(spectrum[0] == spectrum[1])
@@ -104,16 +103,13 @@ def test_compute_spectrum_medfilt(tsig, tsig_sine):
     freqs, spectrum = compute_spectrum_medfilt(tsig, FS)
     assert freqs.shape == spectrum.shape
 
-    # Compute raw estimate of psd using fourier transform
-    #   Only look at the spectrum up to the Nyquist frequency
+    # Compute raw estimate of psd using FFT
     sig_len = len(tsig_sine)
-    nyq_freq = sig_len//2
-    sig_ft = np.fft.fft(tsig_sine)[:nyq_freq]
-    psd = np.abs(sig_ft)**2/(FS * sig_len)
+    psd = np.abs(np.fft.rfft(tsig_sine))**2 / (FS * sig_len)
 
     # The medfilt here should be taking the median of a window with one sample
     #   Therefore, it should match the estimate of psd from above
-    _, psd_medfilt = compute_spectrum(tsig_sine, FS, method='medfilt', filt_len=0.1)
+    _, psd_medfilt = compute_spectrum_medfilt(tsig_sine, FS, filt_len=0.1)
     assert np.allclose(psd, psd_medfilt, atol=EPS)
 
 def test_compute_spectrum_multitaper(tsig_sine, tsig2d):
@@ -131,4 +127,3 @@ def test_compute_spectrum_multitaper(tsig_sine, tsig2d):
     idx_freq_sine = np.argmin(np.abs(freqs - FREQ_SINE))
     idx_peak = np.argmax(spectrum)
     assert idx_freq_sine == idx_peak
-

From f64cd1f7e40bc7dedaf2c2e390aa9371735329b8 Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Sun, 12 Oct 2025 16:39:41 +0100
Subject: [PATCH 4/8] add window option to fft

---
 neurodsp/spectral/power.py            | 12 ++++++++++--
 neurodsp/tests/spectral/test_power.py |  8 ++++++--
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py
index 972992a64..92c865a55 100644
--- a/neurodsp/spectral/power.py
+++ b/neurodsp/spectral/power.py
@@ -9,6 +9,7 @@
 import numpy as np
 from scipy.signal import spectrogram, medfilt
 from scipy.fft import next_fast_len
+from scipy.signal.windows import get_window
 
 from neurodsp.utils.core import get_avg_func
 from neurodsp.utils.data import create_freqs
@@ -72,7 +73,7 @@ def compute_spectrum(sig, fs, method='welch', **kwargs):
 
 SPECTRUM_INPUTS = {
     'wavelet' : ['freqs', 'avg_type', 'n_cycles', 'scaling', 'norm'],
-    'fft' : ['f_range'],
+    'fft' : ['window', 'f_range'],
     'welch' : ['avg_type', 'window', 'nperseg', 'noverlap', 'nfft', \
                'fast_len', 'f_range', 'outlier_percent'],
     'medfilt' : ['filt_len', 'f_range'],
@@ -139,7 +140,7 @@ def compute_spectrum_wavelet(sig, fs, freqs, avg_type='mean', **kwargs):
 
 
 @multidim(select=[0])
-def compute_spectrum_fft(sig, fs, f_range=None):
+def compute_spectrum_fft(sig, fs, window=None, f_range=None):
     """Compute the power spectrum based on a single FFT.
 
     Parameters
@@ -148,6 +149,10 @@ def compute_spectrum_fft(sig, fs, f_range=None):
         Time series.
     fs : float
         Sampling rate, in Hz.
+    window : str, tuple, float
+        Window function to apply to signal.
+        Typically, this is a string of the name of the window to use (e.g. 'hann' or 'hamming').
+        See `scipy.signal.windows.get_window` for details.
     f_range : list of [float, float], optional
         Frequency range to sub-select from the power spectrum.
 
@@ -159,6 +164,9 @@ def compute_spectrum_fft(sig, fs, f_range=None):
         Power spectral density.
     """
 
+    if window is not None:
+        sig = sig * get_window(window, len(sig))
+
     # Compute the FFT and convert to power & compute corresponding frequency vector
     spectrum = np.abs(np.fft.rfft(sig)) ** 2.
     freqs = np.fft.rfftfreq(len(sig), 1. / fs)
diff --git a/neurodsp/tests/spectral/test_power.py b/neurodsp/tests/spectral/test_power.py
index 5df0a3342..791a5fd14 100644
--- a/neurodsp/tests/spectral/test_power.py
+++ b/neurodsp/tests/spectral/test_power.py
@@ -57,8 +57,12 @@ def test_compute_spectrum_2d(tsig2d):
 
 def test_compute_spectrum_fft(tsig, tsig_sine):
 
-    freqs, spectrum = compute_spectrum_fft(tsig, FS)
-    assert freqs.shape == spectrum.shape
+    freqs1, spectrum1 = compute_spectrum_fft(tsig, FS)
+    assert freqs1.shape == spectrum1.shape
+
+    # Test applying a window function
+    freqs2, spectrum2 = compute_spectrum_fft(tsig, FS, window='hann')
+    assert freqs2.shape == spectrum2.shape
 
 def test_compute_spectrum_welch(tsig, tsig_sine):
 

From afdca43cd5da1b7339b852a254bcaf34c2b5cc14 Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Sun, 12 Oct 2025 17:12:09 +0100
Subject: [PATCH 5/8] add & use pad_signal

---
 neurodsp/spectral/power.py            | 12 ++++++++---
 neurodsp/spectral/utils.py            | 30 +++++++++++++++++++++++++++
 neurodsp/tests/spectral/test_power.py |  5 +++++
 neurodsp/tests/spectral/test_utils.py | 22 ++++++++++++++++++++
 4 files changed, 66 insertions(+), 3 deletions(-)

diff --git a/neurodsp/spectral/power.py b/neurodsp/spectral/power.py
index 92c865a55..b9900f958 100644
--- a/neurodsp/spectral/power.py
+++ b/neurodsp/spectral/power.py
@@ -17,7 +17,7 @@
 from neurodsp.utils.checks import check_param_options
 from neurodsp.utils.outliers import discard_outliers
 from neurodsp.timefrequency.wavelets import compute_wavelet_transform
-from neurodsp.spectral.utils import trim_spectrum, window_pad, get_positive_fft_outputs
+from neurodsp.spectral.utils import trim_spectrum, window_pad, pad_signal
 from neurodsp.spectral.checks import check_spg_settings, check_mt_settings
 
 ###################################################################################################
@@ -140,7 +140,7 @@ def compute_spectrum_wavelet(sig, fs, freqs, avg_type='mean', **kwargs):
 
 
 @multidim(select=[0])
-def compute_spectrum_fft(sig, fs, window=None, f_range=None):
+def compute_spectrum_fft(sig, fs, window=None, nfft=None, f_range=None):
     """Compute the power spectrum based on a single FFT.
 
     Parameters
@@ -149,10 +149,13 @@ def compute_spectrum_fft(sig, fs, window=None, f_range=None):
         Time series.
     fs : float
         Sampling rate, in Hz.
-    window : str, tuple, float
+    window : str or tuple or float, optional
         Window function to apply to signal.
         Typically, this is a string of the name of the window to use (e.g. 'hann' or 'hamming').
         See `scipy.signal.windows.get_window` for details.
+    nfft : int, optional
+        Number of samples per for the FFT estimation.
+        If provided and nfft > len(sig), then the signal is zero-padded to this length.
     f_range : list of [float, float], optional
         Frequency range to sub-select from the power spectrum.
 
@@ -167,6 +170,9 @@ def compute_spectrum_fft(sig, fs, window=None, f_range=None):
     if window is not None:
         sig = sig * get_window(window, len(sig))
 
+    if nfft is not None:
+        sig = pad_signal(sig, nfft)
+
     # Compute the FFT and convert to power & compute corresponding frequency vector
     spectrum = np.abs(np.fft.rfft(sig)) ** 2.
     freqs = np.fft.rfftfreq(len(sig), 1. / fs)
diff --git a/neurodsp/spectral/utils.py b/neurodsp/spectral/utils.py
index 051a8e395..9c251c538 100644
--- a/neurodsp/spectral/utils.py
+++ b/neurodsp/spectral/utils.py
@@ -167,6 +167,36 @@ def get_positive_fft_outputs(freqs, powers=None, drop_zero=False):
         return freqs[start_ind:end_ind]
 
 
+def pad_signal(sig, length):
+    """Pad a signal to a desired length.
+
+    Parameters
+    ----------
+    sig : 1d array
+        Signal to pad.
+    length : int
+        Output length to pad the signal to.
+
+    Returns
+    -------
+    sig : 1d array
+        Padded signal.
+
+    Notes
+    -----
+    This approach pads the signal evenly on the left and right side with 0s.
+    If the padding length ends up being odd, this approach will split to padding to
+    have one less at the front / left side pad, and one more on the right / end side pad.
+    """
+
+    if length > len(sig):
+        npad_total = length - len(sig)
+        npad_left, npad_right = int(np.floor(npad_total / 2)), int(np.ceil(npad_total / 2))
+        sig = np.pad(sig, (npad_left, npad_right), mode='constant', constant_values=0)
+
+    return sig
+
+
 def window_pad(sig, nperseg, noverlap, npad, fast_len,
                nwindows=None, nsamples=None, pad_left=None, pad_right=None):
     """Pads windows (for Welch's PSD) with zeros.
diff --git a/neurodsp/tests/spectral/test_power.py b/neurodsp/tests/spectral/test_power.py
index 791a5fd14..5c8dd961a 100644
--- a/neurodsp/tests/spectral/test_power.py
+++ b/neurodsp/tests/spectral/test_power.py
@@ -64,6 +64,11 @@ def test_compute_spectrum_fft(tsig, tsig_sine):
     freqs2, spectrum2 = compute_spectrum_fft(tsig, FS, window='hann')
     assert freqs2.shape == spectrum2.shape
 
+    # Test padding signal
+    freqs3, spectrum3 = compute_spectrum_fft(tsig, FS, nfft=1.5*len(tsig))
+    assert freqs3.shape == spectrum3.shape
+    assert freqs2.shape != freqs3.shape
+
 def test_compute_spectrum_welch(tsig, tsig_sine):
 
     freqs, spectrum = compute_spectrum_welch(tsig, FS, avg_type='mean')
diff --git a/neurodsp/tests/spectral/test_utils.py b/neurodsp/tests/spectral/test_utils.py
index 7a6ec0e10..dd29ce09a 100644
--- a/neurodsp/tests/spectral/test_utils.py
+++ b/neurodsp/tests/spectral/test_utils.py
@@ -65,6 +65,28 @@ def test_get_positive_fft_outputs():
     assert freqs_out_even.shape == powers_out_even.shape
     assert np.array_equal(freqs_out_even, freqs_even[:3])
 
+def test_pad_signal():
+
+    # Test case: odd length, even number added per side
+    length = 5
+    out1 = pad_signal(np.array([1, 2, 3]), length)
+    assert len(out1) == length
+
+    # Test case: even length, even number added per side
+    length = 6
+    out2 = pad_signal(np.array([1, 2]), length)
+    assert len(out2) == length
+
+    # Test case: odd length, uneven number added per side
+    length = 5
+    out3 = pad_signal(np.array([1, 2]), length)
+    assert len(out3) == length
+
+    # Test case: even length, uneven number added per side
+    length = 6
+    out4 = pad_signal(np.array([1, 2, 3]), length)
+    assert len(out4) == length
+
 @pytest.mark.parametrize("fast_len", [True, False])
 def test_window_pad(fast_len):
 

From 7040e906c7c5f043b45ed0ee45a40a97b3064390 Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Mon, 27 Oct 2025 19:38:27 +0000
Subject: [PATCH 6/8] add fast_len to pad_signal

---
 neurodsp/spectral/utils.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/neurodsp/spectral/utils.py b/neurodsp/spectral/utils.py
index b09e0d0a5..f0883bd58 100644
--- a/neurodsp/spectral/utils.py
+++ b/neurodsp/spectral/utils.py
@@ -167,7 +167,7 @@ def get_positive_fft_outputs(freqs, powers=None, drop_zero=False):
         return freqs[start_ind:end_ind]
 
 
-def pad_signal(sig, length):
+def pad_signal(sig, length, fast_len=False):
     """Pad a signal to a desired length.
 
     Parameters
@@ -176,6 +176,9 @@ def pad_signal(sig, length):
         Signal to pad.
     length : int
         Output length to pad the signal to.
+    fast_len : bool, optional, default: False
+        If True, updates length to the next fastest length to reduce computation time.
+        See scipy.fft.next_fast_len for details.
 
     Returns
     -------
@@ -190,6 +193,8 @@ def pad_signal(sig, length):
     """
 
     if length > len(sig):
+        if fast_len:
+            length = next_fast_len(length)
         npad_total = length - len(sig)
         npad_left, npad_right = int(np.floor(npad_total / 2)), int(np.ceil(npad_total / 2))
         sig = np.pad(sig, (npad_left, npad_right), mode='constant', constant_values=0)

From 445487208948f18fe01ca81936dd213e95902f31 Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Wed, 29 Oct 2025 10:57:46 +0000
Subject: [PATCH 7/8] add compute_spectrum_fft to init & api list

---
 doc/api.rst                   | 1 +
 neurodsp/spectral/__init__.py | 5 +++--
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/doc/api.rst b/doc/api.rst
index e76e040bd..4b7464b5a 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -111,6 +111,7 @@ Spectral Power
     :toctree: generated/
 
     compute_spectrum
+    compute_spectrum_fft
     compute_spectrum_welch
     compute_spectrum_wavelet
     compute_spectrum_medfilt
diff --git a/neurodsp/spectral/__init__.py b/neurodsp/spectral/__init__.py
index c63c4e419..74d81b904 100644
--- a/neurodsp/spectral/__init__.py
+++ b/neurodsp/spectral/__init__.py
@@ -1,7 +1,8 @@
 """Spectral module, for calculating power spectra, spectral variance, etc."""
 
-from .power import (compute_spectrum, compute_spectrum_welch, compute_spectrum_wavelet,
-                    compute_spectrum_medfilt, compute_spectrum_multitaper)
+from .power import (compute_spectrum, compute_spectrum_fft, compute_spectrum_welch,
+                    compute_spectrum_wavelet, compute_spectrum_medfilt,
+                    compute_spectrum_multitaper)
 from .measures import compute_absolute_power, compute_relative_power, compute_band_ratio
 from .variance import compute_scv, compute_scv_rs, compute_spectral_hist
 from .utils import trim_spectrum, trim_spectrogram

From 423f1330f6d80e1da286dad450b723d064e117c9 Mon Sep 17 00:00:00 2001
From: Tom Donoghue <tdonoghue.research@gmail.com>
Date: Wed, 29 Oct 2025 11:00:05 +0000
Subject: [PATCH 8/8] drop get_positive_fft_outputs func

---
 neurodsp/spectral/utils.py            | 37 ---------------------------
 neurodsp/tests/spectral/test_utils.py | 22 ----------------
 2 files changed, 59 deletions(-)

diff --git a/neurodsp/spectral/utils.py b/neurodsp/spectral/utils.py
index f0883bd58..5db1454e6 100644
--- a/neurodsp/spectral/utils.py
+++ b/neurodsp/spectral/utils.py
@@ -130,43 +130,6 @@ def trim_spectrogram(freqs, times, spg, f_range=None, t_range=None):
     return freqs_ext, times_ext, spg_ext
 
 
-def get_positive_fft_outputs(freqs, powers=None, drop_zero=False):
-    """Get the positive frequency values for an FFT.
-
-    Parameters
-    ----------
-    freqs : 1d array
-        Frequency vector corresponding to the FFT estimate, with positive & negative frequencies.
-    powers : 1d array, optional
-        Complex power value estimates from the FFT.
-    drop_zero : bool, optional, default: False
-        Whether to drop the estimate for frequency of 0.
-
-    Returns
-    -------
-    freqs : 1d array
-        Frequencies at which the measure was calculated.
-    spectrum : array
-        Power spectral density.
-        Only returned if an input power spectrum is passed.
-
-    Notes
-    -----
-    This can be used to extract positive only frequency from an FFT, for example,
-    as returned by `np.fft.fft` & np.fft.fftfreq`.
-    """
-
-    start_ind = 1 if drop_zero else 0
-
-    # Get the max positive ind as half length, rounded up if the length is odd
-    end_ind = int(np.ceil(len(freqs) / 2))
-
-    if powers is not None:
-        return freqs[start_ind:end_ind], powers[start_ind:end_ind]
-    else:
-        return freqs[start_ind:end_ind]
-
-
 def pad_signal(sig, length, fast_len=False):
     """Pad a signal to a desired length.
 
diff --git a/neurodsp/tests/spectral/test_utils.py b/neurodsp/tests/spectral/test_utils.py
index 8063d6bf6..8585dec6d 100644
--- a/neurodsp/tests/spectral/test_utils.py
+++ b/neurodsp/tests/spectral/test_utils.py
@@ -43,28 +43,6 @@ def test_trim_spectrogram():
     assert_equal(f_ext, np.array([6, 7, 8]))
     assert_equal(t_ext, times)
 
-def test_get_positive_fft_outputs():
-
-    # Test odd length
-    freqs_odd = np.array([0, 1, 2, -2, -1])
-    powers_odd = np.array([0, 1, 2, 3, 4])
-    freqs_out_odd, powers_out_odd = get_positive_fft_outputs(freqs_odd, powers_odd)
-    assert freqs_out_odd.shape == powers_out_odd.shape
-    assert np.array_equal(freqs_out_odd, freqs_odd[:3])
-
-    # Test drop zero
-    freqs_out_odd2, powers_out_odd2 = get_positive_fft_outputs(freqs_odd, powers_odd, drop_zero=True)
-    assert freqs_out_odd2[0] != 0
-    assert powers_out_odd2[0] != 0
-
-    # Test even length
-    freqs_even = np.array([0, 1, 2, -3, -2, -1])
-    powers_even = np.array([0, 1, 2, 3, 4, 5])
-    freqs_out_even, powers_out_even = get_positive_fft_outputs(freqs_odd, powers_odd)
-    assert freqs_out_even.shape == powers_out_even.shape
-    assert freqs_out_even.shape == powers_out_even.shape
-    assert np.array_equal(freqs_out_even, freqs_even[:3])
-
 def test_pad_signal():
 
     # Test case: odd length, even number added per side