Spectral

Spectral analysis functions: FFT amplitude spectrum, power spectral density, cross-spectral density, coherence, autocorrelation, and cross-correlation.

Three estimators of the same quantity live here, and they fail differently. psd averages periodograms and is the default. blackman_tukey_psd transforms a tapered autocorrelation and trades resolution for variance through the lag window, explicitly. ar_psd fits an all-pole model to the whole record and can separate two close peaks on a record too short for either of the others — at the cost of an order that decides what the answer says.

segment_advice answers the question the first of those poses. nperseg can be wrong in two opposite directions, and only one of them is obvious.

All estimators use scipy.signal under the hood with engineering-friendly defaults (Hann window, density scaling, mean detrending).

---

dspkit.spectral.fft_spectrum(x, fs, window='hann', scaling='amplitude')

Single-sided FFT amplitude spectrum with window amplitude correction.

For a pure sine of amplitude A at frequency f, the returned spectrum will show A at that frequency bin (with scaling='amplitude').

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	Time-domain signal.	required
fs	float	Sampling frequency [Hz].	required
window	str or None	Window function name accepted by scipy.signal.get_window. None uses a rectangular window (no windowing).	'hann'
scaling	(amplitude, rms)	'amplitude' returns peak amplitude per bin. 'rms' returns RMS amplitude (peak / sqrt(2)), useful for comparing sinusoidal components with broadband levels.	'amplitude'

Returns:

Name	Type	Description
freqs	(ndarray, shape(N // 2 + 1))	Frequency vector [Hz].
amplitude	(ndarray, shape(N // 2 + 1))	Amplitude spectrum in the same units as x.

Source code in dspkit/spectral.py

def fft_spectrum(
    x: np.ndarray,
    fs: float,
    window: str | None = "hann",
    scaling: Literal["amplitude", "rms"] = "amplitude",
) -> tuple[np.ndarray, np.ndarray]:
    """
    Single-sided FFT amplitude spectrum with window amplitude correction.

    For a pure sine of amplitude A at frequency f, the returned spectrum will
    show A at that frequency bin (with ``scaling='amplitude'``).

    Parameters
    ----------
    x : array_like, shape (N,)
        Time-domain signal.
    fs : float
        Sampling frequency [Hz].
    window : str or None
        Window function name accepted by ``scipy.signal.get_window``.
        ``None`` uses a rectangular window (no windowing).
    scaling : {'amplitude', 'rms'}
        ``'amplitude'`` returns peak amplitude per bin.
        ``'rms'`` returns RMS amplitude (peak / sqrt(2)), useful for
        comparing sinusoidal components with broadband levels.

    Returns
    -------
    freqs : ndarray, shape (N//2 + 1,)
        Frequency vector [Hz].
    amplitude : ndarray, shape (N//2 + 1,)
        Amplitude spectrum in the same units as ``x``.
    """
    x = np.asarray(x, dtype=float)
    N = len(x)

    if window is not None:
        win = _signal.get_window(window, N)
        # Amplitude correction: scale so that a sine's peak is preserved
        acf = N / win.sum()
        x = x * win
    else:
        acf = 1.0

    X = np.fft.rfft(x)
    freqs = np.fft.rfftfreq(N, d=1.0 / fs)

    amplitude = np.abs(X) * acf / N
    # Double interior bins to account for the discarded negative-frequency mirror
    amplitude[1:-1] *= 2.0

    if scaling == "rms":
        amplitude /= np.sqrt(2.0)

    return freqs, amplitude

---

dspkit.spectral.psd(x, fs, window='hann', nperseg=None, noverlap=None, scaling='density', detrend='constant')

Power spectral density (or power spectrum) via Welch's method.

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	Time-domain signal.	required
fs	float	Sampling frequency [Hz].	required
window	str	Window function (default 'hann').	'hann'
nperseg	int or None	Segment length. Defaults to min(len(x), 1024).	None
noverlap	int or None	Number of overlapping samples between segments. Defaults to nperseg // 2 (50 % overlap).	None
scaling	(density, spectrum)	'density' → PSD [units²/Hz]. 'spectrum' → power spectrum [units²].	'density'
detrend	str or False	Detrending applied to each segment before windowing. 'constant' removes the mean, 'linear' removes a linear trend, False skips detrending.	'constant'

Returns:

Name	Type	Description
freqs	ndarray	Frequency vector [Hz].
Pxx	ndarray	One-sided PSD or power spectrum (real, non-negative).

Source code in dspkit/spectral.py

def psd(
    x: np.ndarray,
    fs: float,
    window: str = "hann",
    nperseg: int | None = None,
    noverlap: int | None = None,
    scaling: Literal["density", "spectrum"] = "density",
    detrend: str | Literal[False] = "constant",
) -> tuple[np.ndarray, np.ndarray]:
    """
    Power spectral density (or power spectrum) via Welch's method.

    Parameters
    ----------
    x : array_like, shape (N,)
        Time-domain signal.
    fs : float
        Sampling frequency [Hz].
    window : str
        Window function (default ``'hann'``).
    nperseg : int or None
        Segment length. Defaults to ``min(len(x), 1024)``.
    noverlap : int or None
        Number of overlapping samples between segments.
        Defaults to ``nperseg // 2`` (50 % overlap).
    scaling : {'density', 'spectrum'}
        ``'density'`` → PSD [units²/Hz].
        ``'spectrum'`` → power spectrum [units²].
    detrend : str or False
        Detrending applied to each segment before windowing.
        ``'constant'`` removes the mean, ``'linear'`` removes a linear trend,
        ``False`` skips detrending.

    Returns
    -------
    freqs : ndarray
        Frequency vector [Hz].
    Pxx : ndarray
        One-sided PSD or power spectrum (real, non-negative).
    """
    x = np.asarray(x, dtype=float)
    if nperseg is None:
        nperseg = min(len(x), 1024)

    freqs, Pxx = _signal.welch(
        x,
        fs=fs,
        window=window,
        nperseg=nperseg,
        noverlap=noverlap,
        scaling=scaling,
        detrend=detrend,
    )
    return freqs, Pxx

---

dspkit.spectral.ar_psd(x, fs, order=None, method='burg', n_freqs=2048, max_order=None, criterion='aic', detrend=True)

Power spectral density from an autoregressive model.

The third classical spectral estimator, alongside Welch (psd) and the correlogram (blackman_tukey_psd). Instead of averaging periodograms it fits an all-pole model to the record and evaluates its transfer function::

S(f) = e / | 1 + sum_k a_k exp(-2 pi i f k / fs) |^2 / fs

Reach for it when the record is too short for Welch. Welch's resolution is set by the segment length and its variance by the number of segments, and a short record cannot give both. An AR model is not segmented at all: it spends the whole record on one fit, so it can resolve two close modes where Welch would need more data than exists.

Measured on two tones at 10.0 and 10.8 Hz in noise, 2 s at 100 Hz -- 200 samples in total::

Welch, nperseg=200 (Be = 0.75 Hz)   one peak at 10.50 Hz
Welch, nperseg=128 (Be = 1.17 Hz)   one peak at 10.16 Hz
Welch, nperseg= 64 (Be = 2.34 Hz)   one peak at 10.94 Hz
AR Burg, order 30                   9.99 and 10.85 Hz

Welch merges them at every usable segment length; the AR fit separates them to within 0.05 Hz. Note the order needed, which is the next point.

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	The record.	required
fs	float	Sampling frequency [Hz].	required
order	int or None	Model order. None selects one by criterion up to max_order.	None
method	(burg, yule_walker)	Burg by default -- it does not window the data and always returns a stable model. Yule-Walker is offered because it is what most textbooks derive and is useful for comparison.	'burg'
n_freqs	int	Frequencies at which to evaluate, from 0 to Nyquist. This is a drawing resolution, not an information one: the model has order poles however finely it is sampled.	2048
max_order	int or None	Ceiling for automatic selection. Defaults to min(N // 2, 100).	None
criterion	(aic, bic)	Used only when order is None.	'aic'
detrend	bool	Remove the mean first (default True). A DC offset is a pole at zero frequency and will otherwise dominate the fit.	True

Returns:

Name	Type	Description
freqs	ndarray	Frequency vector [Hz].
pxx	ndarray	One-sided PSD [units^2/Hz].
info	dict	order, method, criterion, variance (the residual prediction error), reflection_stable and, when the order was chosen, scores.

Notes

The order is the estimate. Unlike Welch, where a bad nperseg gives a blurred but honest picture, a bad AR order changes what the spectrum says. Too low and close modes merge into one broad peak; too high and the model spends poles on noise. Measured on the library's 2-DOF example, displacement output (N = 20480, fs = 1024 Hz, true modes 8.613 Hz at 1.22% damping and 20.795 Hz at 2.94%)::

order   2    one peak, at 10.5 Hz — the two modes merged
order   4    one peak, at 10.5 Hz
order  22    8.504 and 21.761 Hz
order 100    8.504 and 20.760 Hz
AIC -> 83    8.504 and 20.760 Hz
BIC -> 46    8.504 and 21.010 Hz

An order of 4 is not enough for a 2-DOF system, which is the trap: the rule of thumb "two poles per mode" describes a noiseless model, not a fit to a finite noisy record. AIC and BIC both landed somewhere workable here, but AIC will happily run to whatever max_order allows on a long record — it chose the ceiling on the acceleration channel of the same fixture.

It is a model, not a measurement. An AR spectrum is smooth and confident everywhere, including where the data says nothing, and it has no equivalent of a confidence interval falling out of the segment count. Cross-check against psd before believing a peak that only the AR estimate shows.

Peaks are also biased: an all-pole model represents a resonance exactly and an antiresonance only by cancellation, so notches come out shallower than they are. Use psd where the notch matters.

See Also

psd, blackman_tukey_psd, segment_advice, ar_order_selection

Examples:

>>> import numpy as np
>>> from dspkit.spectral import ar_psd
>>> from dspkit._testing import generate_2dof
>>> t, a1, a2 = generate_2dof(duration=20.0, fs=1024.0, seed=0)
>>> f, p, info = ar_psd(a1, 1024.0, order=40)
>>> info['order'], info['reflection_stable']
(40, True)

Source code in dspkit/spectral.py

def ar_psd(
    x: np.ndarray,
    fs: float,
    order: int | None = None,
    method: str = "burg",
    n_freqs: int = 2048,
    max_order: int | None = None,
    criterion: str = "aic",
    detrend: bool = True,
) -> tuple[np.ndarray, np.ndarray, dict]:
    """
    Power spectral density from an autoregressive model.

    The third classical spectral estimator, alongside Welch (`psd`) and the
    correlogram (`blackman_tukey_psd`). Instead of averaging periodograms it
    fits an all-pole model to the record and evaluates its transfer function::

        S(f) = e / | 1 + sum_k a_k exp(-2 pi i f k / fs) |^2 / fs

    **Reach for it when the record is too short for Welch.** Welch's resolution
    is set by the segment length and its variance by the number of segments,
    and a short record cannot give both. An AR model is not segmented at all:
    it spends the whole record on one fit, so it can resolve two close modes
    where Welch would need more data than exists.

    Measured on two tones at 10.0 and 10.8 Hz in noise, 2 s at 100 Hz -- 200
    samples in total::

        Welch, nperseg=200 (Be = 0.75 Hz)   one peak at 10.50 Hz
        Welch, nperseg=128 (Be = 1.17 Hz)   one peak at 10.16 Hz
        Welch, nperseg= 64 (Be = 2.34 Hz)   one peak at 10.94 Hz
        AR Burg, order 30                   9.99 and 10.85 Hz

    Welch merges them at every usable segment length; the AR fit separates them
    to within 0.05 Hz. Note the order needed, which is the next point.

    Parameters
    ----------
    x : array_like, shape (N,)
        The record.
    fs : float
        Sampling frequency [Hz].
    order : int or None
        Model order. ``None`` selects one by `criterion` up to `max_order`.
    method : {'burg', 'yule_walker'}
        Burg by default -- it does not window the data and always returns a
        stable model. Yule-Walker is offered because it is what most textbooks
        derive and is useful for comparison.
    n_freqs : int
        Frequencies at which to evaluate, from 0 to Nyquist. This is a drawing
        resolution, not an information one: the model has `order` poles however
        finely it is sampled.
    max_order : int or None
        Ceiling for automatic selection. Defaults to ``min(N // 2, 100)``.
    criterion : {'aic', 'bic'}
        Used only when `order` is None.
    detrend : bool
        Remove the mean first (default True). A DC offset is a pole at zero
        frequency and will otherwise dominate the fit.

    Returns
    -------
    freqs : ndarray
        Frequency vector [Hz].
    pxx : ndarray
        One-sided PSD [units^2/Hz].
    info : dict
        ``order``, ``method``, ``criterion``, ``variance`` (the residual
        prediction error), ``reflection_stable`` and, when the order was
        chosen, ``scores``.

    Notes
    -----
    **The order is the estimate.** Unlike Welch, where a bad `nperseg` gives a
    blurred but honest picture, a bad AR order changes what the spectrum *says*.
    Too low and close modes merge into one broad peak; too high and the model
    spends poles on noise. Measured on the library's 2-DOF example, displacement
    output (N = 20480, fs = 1024 Hz, true modes 8.613 Hz at 1.22% damping and
    20.795 Hz at 2.94%)::

        order   2    one peak, at 10.5 Hz — the two modes merged
        order   4    one peak, at 10.5 Hz
        order  22    8.504 and 21.761 Hz
        order 100    8.504 and 20.760 Hz
        AIC -> 83    8.504 and 20.760 Hz
        BIC -> 46    8.504 and 21.010 Hz

    An order of 4 is *not* enough for a 2-DOF system, which is the trap: the
    rule of thumb "two poles per mode" describes a noiseless model, not a fit to
    a finite noisy record. AIC and BIC both landed somewhere workable here, but
    AIC will happily run to whatever `max_order` allows on a long record — it
    chose the ceiling on the acceleration channel of the same fixture.

    **It is a model, not a measurement.** An AR spectrum is smooth and confident
    everywhere, including where the data says nothing, and it has no equivalent
    of a confidence interval falling out of the segment count. Cross-check
    against `psd` before believing a peak that only the AR estimate shows.

    Peaks are also biased: an all-pole model represents a resonance exactly and
    an antiresonance only by cancellation, so notches come out shallower than
    they are. Use `psd` where the notch matters.

    See Also
    --------
    psd, blackman_tukey_psd, segment_advice, ar_order_selection

    Examples
    --------
    >>> import numpy as np
    >>> from dspkit.spectral import ar_psd
    >>> from dspkit._testing import generate_2dof
    >>> t, a1, a2 = generate_2dof(duration=20.0, fs=1024.0, seed=0)
    >>> f, p, info = ar_psd(a1, 1024.0, order=40)
    >>> info['order'], info['reflection_stable']
    (40, True)
    """
    if method not in AR_METHODS:
        raise ValueError(
            f"method must be one of {AR_METHODS}, not {method!r}."
        )
    x = np.asarray(x, dtype=float).ravel()
    if x.size < 8:
        raise ValueError(f"An AR fit needs more than {x.size} samples.")
    if detrend:
        x = x - x.mean()

    scores = None
    if order is None:
        if max_order is None:
            max_order = int(min(x.size // 2, 100))
        order, scores = ar_order_selection(x, max_order, method, criterion)
    order = int(order)
    if order < 1:
        raise ValueError("order must be at least 1.")
    if order >= x.size:
        raise ValueError(
            f"order={order} needs more than {x.size} samples to fit."
        )

    a, e = (_burg if method == "burg" else _yule_walker)(x, order)

    freqs = np.linspace(0.0, fs / 2.0, int(n_freqs))
    k = np.arange(order + 1)[:, None]
    denom = a[:, None] * np.exp(-2j * np.pi * freqs[None, :] * k / fs)
    h = np.abs(denom.sum(axis=0)) ** 2
    # One-sided: the two-sided density e/(fs |A|^2) doubled everywhere except
    # DC and Nyquist, which have no mirror image to fold in.
    pxx = e / (fs * np.maximum(h, np.finfo(float).tiny))
    pxx[1:-1] *= 2.0

    roots = np.roots(a) if order > 0 else np.array([])
    info = {
        "order": order,
        "method": method,
        "criterion": criterion if scores is not None else None,
        "variance": float(e),
        "reflection_stable": bool(np.all(np.abs(roots) < 1.0)) if roots.size else True,
    }
    if scores is not None:
        info["scores"] = scores
    return freqs, pxx, info

---

dspkit.spectral.ar_order_selection(x, max_order, method='burg', criterion='aic')

Choose an AR order by an information criterion.

Returns (best_order, scores) where scores[p] is the criterion at order p. "aic" is the Akaike criterion N ln(e_p) + 2p, "bic" swaps the penalty for p ln N and so picks lower orders.

Order selection on a spectrum is not the same problem as on a forecast, and no criterion settles it: too low and peaks merge, too high and the estimate grows spurious ones. Treat the answer as a starting point and look at the spectrum.

Source code in dspkit/spectral.py

def ar_order_selection(
    x: np.ndarray,
    max_order: int,
    method: str = "burg",
    criterion: str = "aic",
) -> tuple[int, np.ndarray]:
    """
    Choose an AR order by an information criterion.

    Returns ``(best_order, scores)`` where ``scores[p]`` is the criterion at
    order ``p``. ``"aic"`` is the Akaike criterion ``N ln(e_p) + 2p``,
    ``"bic"`` swaps the penalty for ``p ln N`` and so picks lower orders.

    Order selection on a *spectrum* is not the same problem as on a forecast,
    and no criterion settles it: too low and peaks merge, too high and the
    estimate grows spurious ones. Treat the answer as a starting point and look
    at the spectrum.
    """
    x = np.asarray(x, dtype=float).ravel()
    n = x.size
    fit = _burg if method == "burg" else _yule_walker
    # One pass. Both recursions produce the order-p error on their way to
    # order `max_order`, so refitting per candidate did the same arithmetic
    # `max_order` times over -- measurable on any real record, and the reason
    # the app's AR tab used to sit there computing.
    _, _, errs = fit(x, int(max_order), history=True)
    p_idx = np.arange(max_order + 1, dtype=float)
    with np.errstate(divide="ignore", invalid="ignore"):
        log_e = np.log(np.where(errs > 0, errs, np.nan))
    penalty = 2.0 * p_idx if criterion == "aic" else p_idx * np.log(n)
    scores = n * log_e + penalty
    scores[0] = np.inf                      # order 0 is not a model
    scores = np.where(np.isfinite(scores), scores, np.inf)
    best = int(np.argmin(scores))
    return max(best, 1), scores

---

dspkit.spectral.csd(x, y, fs, window='hann', nperseg=None, noverlap=None, detrend='constant')

Cross-spectral density via Welch's method.

Gxy(f) = E[X*(f) Y(f)] / Hz, where X, Y are the DFTs of x and y.

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	Input signals. They must have the same sampling frequency.	required
y	(array_like, shape(N))	Input signals. They must have the same sampling frequency.	required
fs	float	Sampling frequency [Hz].	required
window	str	Window function (default 'hann').	'hann'
nperseg	int or None	Segment length. Defaults to min(len(x), len(y), 1024).	None
noverlap	int or None	Overlapping samples. Defaults to nperseg // 2.	None
detrend	str or False	Per-segment detrending (see psd).	'constant'

Returns:

Name	Type	Description
freqs	ndarray	Frequency vector [Hz].
Pxy	ndarray(complex)	One-sided cross-spectral density [units_x · units_y / Hz].

Source code in dspkit/spectral.py

def csd(
    x: np.ndarray,
    y: np.ndarray,
    fs: float,
    window: str = "hann",
    nperseg: int | None = None,
    noverlap: int | None = None,
    detrend: str | Literal[False] = "constant",
) -> tuple[np.ndarray, np.ndarray]:
    """
    Cross-spectral density via Welch's method.

    Gxy(f) = E[X*(f) Y(f)] / Hz,  where X, Y are the DFTs of x and y.

    Parameters
    ----------
    x, y : array_like, shape (N,)
        Input signals. They must have the same sampling frequency.
    fs : float
        Sampling frequency [Hz].
    window : str
        Window function (default ``'hann'``).
    nperseg : int or None
        Segment length. Defaults to ``min(len(x), len(y), 1024)``.
    noverlap : int or None
        Overlapping samples. Defaults to ``nperseg // 2``.
    detrend : str or False
        Per-segment detrending (see `psd`).

    Returns
    -------
    freqs : ndarray
        Frequency vector [Hz].
    Pxy : ndarray (complex)
        One-sided cross-spectral density [units_x · units_y / Hz].
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if nperseg is None:
        nperseg = min(len(x), len(y), 1024)

    freqs, Pxy = _signal.csd(
        x,
        y,
        fs=fs,
        window=window,
        nperseg=nperseg,
        noverlap=noverlap,
        detrend=detrend,
    )
    return freqs, Pxy

---

dspkit.spectral.coherence(x, y, fs, window='hann', nperseg=None, noverlap=None, detrend='constant', min_segments=8)

Magnitude-squared coherence between x and y.

Cxy(f) = |Gxy(f)|² / (Gxx(f) · Gyy(f)), values in [0, 1].

A value near 1 means the two signals are linearly related at that frequency. A value near 0 indicates noise or nonlinearity.

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))		required
y	(array_like, shape(N))		required
fs	float	Sampling frequency [Hz].	required
window	str	Window function (default 'hann').	'hann'
nperseg	int or None	Segment length. Defaults to min(len(x), len(y), 1024).	None
noverlap	int or None	Overlapping samples. Defaults to nperseg // 2.	None
detrend	str or False	Per-segment detrending.	'constant'
min_segments	int	Warn below this many Welch segments (default 8). Fewer than two segments raises instead — see Notes.	8

Returns:

Name	Type	Description
freqs	ndarray
Cxy	ndarray	Magnitude-squared coherence, values in [0, 1].

Raises:

Type	Description
ValueError	If the parameters give fewer than two Welch segments.

Notes

Coherence is made by the averaging, not by the formula. Within a single segment |Gxy|² = Gxx·Gyy identically, so a one-segment estimate is exactly 1.0 at every frequency whatever the two signals are. nperseg = len(x) is precisely that case, and it is rejected rather than returned. Measured on the library's 2-DOF example (N = 20480, fs = 1024 Hz): mean coherence 0.0675 at nperseg=1024 (39 segments), and exactly 1.0000 everywhere at nperseg=N.

The bias does not stop at one segment, it only becomes finite. For two independent signals averaged over n_d segments the expected coherence is about 1 / n_d, so a value of 0.2 from 5 segments is what independence looks like, not evidence of a relationship. min_segments warns while that floor is still large; it does not correct for it.

What this will not tell you: whether the relationship is causal, which channel leads (use the cross-spectrum phase), whether a low value means noise or nonlinearity, or whether a high value at one frequency survives conditioning on the other channels in an array — for that see dspkit.multisensor.partial_coherence.

Source code in dspkit/spectral.py

def coherence(
    x: np.ndarray,
    y: np.ndarray,
    fs: float,
    window: str = "hann",
    nperseg: int | None = None,
    noverlap: int | None = None,
    detrend: str | Literal[False] = "constant",
    min_segments: int = 8,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Magnitude-squared coherence between x and y.

    Cxy(f) = |Gxy(f)|² / (Gxx(f) · Gyy(f)),  values in [0, 1].

    A value near 1 means the two signals are linearly related at that frequency.
    A value near 0 indicates noise or nonlinearity.

    Parameters
    ----------
    x, y : array_like, shape (N,)
    fs : float
        Sampling frequency [Hz].
    window : str
        Window function (default ``'hann'``).
    nperseg : int or None
        Segment length. Defaults to ``min(len(x), len(y), 1024)``.
    noverlap : int or None
        Overlapping samples. Defaults to ``nperseg // 2``.
    detrend : str or False
        Per-segment detrending.
    min_segments : int
        Warn below this many Welch segments (default 8). Fewer than two
        segments raises instead — see Notes.

    Returns
    -------
    freqs : ndarray
    Cxy : ndarray
        Magnitude-squared coherence, values in [0, 1].

    Raises
    ------
    ValueError
        If the parameters give fewer than two Welch segments.

    Notes
    -----
    **Coherence is made by the averaging, not by the formula.** Within a
    single segment |Gxy|² = Gxx·Gyy identically, so a one-segment estimate is
    exactly 1.0 at every frequency whatever the two signals are.
    ``nperseg = len(x)`` is precisely that case, and it is rejected rather
    than returned. Measured on the library's 2-DOF example (N = 20480,
    fs = 1024 Hz): mean coherence 0.0675 at ``nperseg=1024`` (39 segments),
    and exactly 1.0000 everywhere at ``nperseg=N``.

    The bias does not stop at one segment, it only becomes finite. For two
    independent signals averaged over ``n_d`` segments the expected coherence
    is about ``1 / n_d``, so a value of 0.2 from 5 segments is what
    independence looks like, not evidence of a relationship. ``min_segments``
    warns while that floor is still large; it does not correct for it.

    What this will not tell you: whether the relationship is causal, which
    channel leads (use the cross-spectrum phase), whether a low value means
    noise or nonlinearity, or whether a high value at one frequency survives
    conditioning on the other channels in an array — for that see
    ``dspkit.multisensor.partial_coherence``.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if nperseg is None:
        nperseg = min(len(x), len(y), 1024)

    _check_welch_segments(
        "coherence",
        min(len(x), len(y)),
        nperseg,
        noverlap,
        hard_min=2,
        min_segments=min_segments,
        hard_reason=(
            "A single-segment estimate is identically 1.0 at every frequency "
            "by construction and says nothing about the signals."
        ),
    )

    freqs, Cxy = _signal.coherence(
        x,
        y,
        fs=fs,
        window=window,
        nperseg=nperseg,
        noverlap=noverlap,
        detrend=detrend,
    )
    return freqs, Cxy

---

dspkit.spectral.resolution_bandwidth(fs, nperseg, window='hann')

Effective resolution bandwidth of a Welch estimate [Hz].

The equivalent noise bandwidth of the window, N sum(w^2) / (sum w)^2 bins, converted to Hz. It is the width of the rectangle that would pass the same noise power, and it is what decides whether a sharp peak is resolved or smeared -- not the bin spacing fs / nperseg, which is finer and flatters the estimate.

Measured: Hann gives exactly 1.5 bins, Hamming 1.363, rectangular 1.0, Blackman 1.727 and flat-top 3.770, independent of length.

Parameters:

Name	Type	Description	Default
fs	float	Sampling frequency [Hz].	required
nperseg	int	Welch segment length in samples.	required
window	str	Window name, as passed to the estimators.	'hann'

Source code in dspkit/spectral.py

def resolution_bandwidth(fs: float, nperseg: int, window: str = "hann") -> float:
    """
    Effective resolution bandwidth of a Welch estimate [Hz].

    The equivalent noise bandwidth of the window, ``N sum(w^2) / (sum w)^2``
    bins, converted to Hz. It is the width of the rectangle that would pass the
    same noise power, and it is what decides whether a sharp peak is resolved
    or smeared -- not the bin spacing ``fs / nperseg``, which is finer and
    flatters the estimate.

    Measured: Hann gives exactly 1.5 bins, Hamming 1.363, rectangular 1.0,
    Blackman 1.727 and flat-top 3.770, independent of length.

    Parameters
    ----------
    fs : float
        Sampling frequency [Hz].
    nperseg : int
        Welch segment length in samples.
    window : str
        Window name, as passed to the estimators.
    """
    w = _signal.get_window(window, int(nperseg))
    enbw_bins = len(w) * np.sum(w ** 2) / (np.sum(w) ** 2)
    return float(enbw_bins * fs / nperseg)

---

dspkit.spectral.segment_advice(x, fs, nperseg=None, noverlap=None, window='hann', n_inputs=1, target_segments=20, peak_prominence_db=6.0, dynamic_range_db=20.0, max_peaks=8)

Whether nperseg is long enough to resolve the peaks and short enough to average.

The two failure modes pull in opposite directions and only one of them is obvious. Too few averages is the familiar one: coherence is biased up by about q / n_d and spectra are noisy, and the cure is a shorter segment. The other is that a segment too short to resolve the narrowest peak lets leakage smear peaks and notches, which biases coherence down and manufactures a residual where none exists -- an error that grows as the square of the resolution bandwidth and is worst exactly where the signal is strongest. Shortening nperseg to win averages walks straight into it.

This returns both numbers together so the trade can be made deliberately rather than one side at a time.

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	The record. Used to find the sharpest peak; pass the channel whose resonances matter.	required
fs	float	Sampling frequency [Hz].	required
nperseg	int or None	Segment length to assess. Defaults to min(N, 1024), the library's own default, which is what makes this useful before choosing one.	None
noverlap	int or None	Overlap in samples, default nperseg // 2.	None
window	str	Window name.	'hann'
n_inputs	int	Number of predictors, for the coherence bias floor q / n_d. 1 for an ordinary spectrum or a single-input coherence.	1
target_segments	int	Averages considered adequate, used only to size the record this would need. 20 puts the bias floor at 0.05.	20
peak_prominence_db	float	How far a peak must stand above its surroundings, in dB, to count as a resonance worth resolving (default 6, a factor of four in power). Lower it on a record whose modes barely clear the noise.	6.0
dynamic_range_db	float	Ignore peaks more than this far below the tallest (default 20 dB). Widens the net on a flat spectrum; narrow it on one with a huge dynamic range where only the top mode is of interest.	20.0
max_peaks	int	Consider only the tallest this many peaks. "The sharpest peak" means the sharpest one that matters, not the narrowest wiggle on the floor.	8

Returns:

Type	Description
dict with keys	nperseg, n_segments, bias_floor, resolution_bw [Hz], peak_freq and peak_bandwidth [Hz] of the sharpest peak found (None if none was), ratio -- resolution_bw / peak_bandwidth, the quantity the classical rule bounds at 1/4 -- resolved, nperseg_for_resolution, duration_for_both [s], verdict (one of "ok", "unresolved", "marginal", "too_few", "squeezed") and advice, a sentence naming what to change.

Notes

An unresolved peak cannot report its own width. The bandwidth measured off a smeared peak is roughly the resolution bandwidth itself, whatever the true width is, so ratio saturates near 1 rather than growing. That is why the verdict for ratio above about 0.9 is "unresolved" rather than a number: the honest answer is that the record cannot yet say how narrow the peak is, and the test is to lengthen nperseg until the measured width stops shrinking.

"squeezed" means both cannot be satisfied on this record: resolving the peak leaves too few segments. That is a statement about the record, not about the parameters, and the fix is more data. duration_for_both says how much.

See Also

resolution_bandwidth, coherence, dspkit.frf.error_spectrum, dspkit.peaks.peak_bandwidth

Source code in dspkit/spectral.py

def segment_advice(
    x: np.ndarray,
    fs: float,
    nperseg: int | None = None,
    noverlap: int | None = None,
    window: str = "hann",
    n_inputs: int = 1,
    target_segments: int = 20,
    peak_prominence_db: float = 6.0,
    dynamic_range_db: float = 20.0,
    max_peaks: int = 8,
) -> dict:
    """
    Whether ``nperseg`` is long enough to resolve the peaks and short enough to
    average.

    **The two failure modes pull in opposite directions and only one of them is
    obvious.** Too few averages is the familiar one: coherence is biased up by
    about ``q / n_d`` and spectra are noisy, and the cure is a shorter segment.
    The other is that a segment too short to resolve the narrowest peak lets
    leakage smear peaks and notches, which biases coherence *down* and
    manufactures a residual where none exists -- an error that grows as the
    square of the resolution bandwidth and is worst exactly where the signal is
    strongest. Shortening ``nperseg`` to win averages walks straight into it.

    This returns both numbers together so the trade can be made deliberately
    rather than one side at a time.

    Parameters
    ----------
    x : array_like, shape (N,)
        The record. Used to find the sharpest peak; pass the channel whose
        resonances matter.
    fs : float
        Sampling frequency [Hz].
    nperseg : int or None
        Segment length to assess. Defaults to ``min(N, 1024)``, the library's
        own default, which is what makes this useful before choosing one.
    noverlap : int or None
        Overlap in samples, default ``nperseg // 2``.
    window : str
        Window name.
    n_inputs : int
        Number of predictors, for the coherence bias floor ``q / n_d``. 1 for
        an ordinary spectrum or a single-input coherence.
    target_segments : int
        Averages considered adequate, used only to size the record this would
        need. 20 puts the bias floor at 0.05.
    peak_prominence_db : float
        How far a peak must stand above its surroundings, in dB, to count as a
        resonance worth resolving (default 6, a factor of four in power).
        Lower it on a record whose modes barely clear the noise.
    dynamic_range_db : float
        Ignore peaks more than this far below the tallest (default 20 dB).
        Widens the net on a flat spectrum; narrow it on one with a huge
        dynamic range where only the top mode is of interest.
    max_peaks : int
        Consider only the tallest this many peaks. "The sharpest peak" means
        the sharpest one that matters, not the narrowest wiggle on the floor.

    Returns
    -------
    dict with keys
        ``nperseg``, ``n_segments``, ``bias_floor``, ``resolution_bw`` [Hz],
        ``peak_freq`` and ``peak_bandwidth`` [Hz] of the sharpest peak found
        (``None`` if none was), ``ratio`` -- ``resolution_bw / peak_bandwidth``,
        the quantity the classical rule bounds at 1/4 -- ``resolved``,
        ``nperseg_for_resolution``, ``duration_for_both`` [s], ``verdict``
        (one of ``"ok"``, ``"unresolved"``, ``"marginal"``, ``"too_few"``,
        ``"squeezed"``) and ``advice``, a sentence naming what to change.

    Notes
    -----
    **An unresolved peak cannot report its own width.** The bandwidth measured
    off a smeared peak is roughly the resolution bandwidth itself, whatever the
    true width is, so ``ratio`` saturates near 1 rather than growing. That is
    why the verdict for ``ratio`` above about 0.9 is "unresolved" rather than a
    number: the honest answer is that the record cannot yet say how narrow the
    peak is, and the test is to lengthen ``nperseg`` until the measured width
    stops shrinking.

    ``"squeezed"`` means both cannot be satisfied on this record: resolving the
    peak leaves too few segments. That is a statement about the *record*, not
    about the parameters, and the fix is more data. ``duration_for_both`` says
    how much.

    See Also
    --------
    resolution_bandwidth, coherence, dspkit.frf.error_spectrum,
    dspkit.peaks.peak_bandwidth
    """
    from .peaks import find_peaks as _find_peaks, peak_bandwidth as _peak_bw

    x = np.asarray(x, dtype=float).ravel()
    n = x.size
    if nperseg is None:
        nperseg = min(n, 1024)
    nperseg = int(min(nperseg, n))
    step = nperseg // 2 if noverlap is None else nperseg - int(noverlap)
    n_seg = _welch_segment_count(n, nperseg, noverlap)

    be = resolution_bandwidth(fs, nperseg, window)
    floor = n_inputs / n_seg if n_seg > 0 else float("inf")

    # The sharpest peak, measured at this very resolution -- which is the point:
    # if it comes back at the resolution bandwidth, it is not resolved.
    #
    # Detected on the dB spectrum with a prominence threshold, because on a
    # linear PSD every noise wiggle is a peak and the *narrowest* of those is
    # always about one resolution bandwidth wide. Taking it would report
    # "unresolved" on any record whatsoever, which is a diagnostic that never
    # says anything. A resonance worth resolving stands well clear of the floor,
    # so prominence in dB is the right filter, and only the tallest few are
    # considered -- "the sharpest peak" means the sharpest one that matters.
    peak_f = peak_bw = None
    try:
        freqs, pxx = psd(x, fs, window=window, nperseg=nperseg, noverlap=noverlap)
        db = 10.0 * np.log10(np.maximum(pxx, np.finfo(float).tiny))
        pf, pv, _ = _find_peaks(freqs, db, prominence=peak_prominence_db,
                                max_peaks=max_peaks)
        # And within `dynamic_range_db` of the tallest. With few averages a PSD
        # is so noisy that 6 dB of prominence is nothing -- measured on a 2-DOF
        # record at 4 averages, the narrowest "significant" peak was a noise
        # spike at 182 Hz, fifty-odd dB below the modes. A resonance that
        # matters sits near the top of the spectrum, and saying so is what keeps
        # this from tracking the noise floor as the resolution improves.
        if len(pf):
            keep = pv >= (np.max(pv) - dynamic_range_db)
            pf, pv = pf[keep], pv[keep]
        if len(pf):
            # Width is measured on the linear PSD: half-power is a factor of two
            # in power, which is not what half-height means in dB.
            pfs, bws, _ = _peak_bw(freqs, pxx, pf)
            good = np.isfinite(bws) & (bws > 0)
            if np.any(good):
                i = int(np.argmin(bws[good]))
                peak_f = float(pfs[good][i])
                peak_bw = float(bws[good][i])
    except (ValueError, IndexError):
        pass

    ratio = (be / peak_bw) if peak_bw else None
    resolved = bool(ratio is not None and ratio <= 0.25)

    # What length would resolve it, and what record length would then still
    # leave enough averages at the caller's overlap fraction.
    nperseg_res = None
    duration_both = None
    if peak_bw:
        enbw_bins = be * nperseg / fs
        nperseg_res = int(np.ceil(4.0 * enbw_bins * fs / peak_bw))
        frac = 1.0 - step / nperseg
        samples = nperseg_res * (target_segments * (1.0 - frac) + frac)
        duration_both = float(samples / fs)

    if peak_bw is None:
        verdict = "ok" if n_seg >= target_segments else "too_few"
    elif ratio > 0.9:
        verdict = "unresolved"
    elif ratio > 0.25:
        verdict = "marginal"
    elif n_seg < target_segments:
        verdict = "too_few"
    else:
        verdict = "ok"

    if verdict in ("unresolved", "marginal") and duration_both and \
            duration_both > n / fs * 1.05:
        verdict = "squeezed"

    if verdict == "ok":
        advice = (
            f"nperseg={nperseg} resolves the sharpest peak "
            f"(Be/Br = {ratio:.2f}) and leaves {n_seg} segments "
            f"(bias floor {floor:.3f}). Nothing to change."
            if peak_bw else
            f"No clear peak to resolve; {n_seg} segments is adequate "
            f"(bias floor {floor:.3f})."
        )
    elif verdict == "unresolved":
        advice = (
            f"The sharpest peak near {peak_f:.4g} Hz is not resolved: its "
            f"apparent width {peak_bw:.4g} Hz is the resolution bandwidth "
            f"{be:.4g} Hz, so its true width could be anything smaller. "
            f"Raise nperseg to at least {nperseg_res} and see whether the width "
            f"keeps shrinking. Do not shorten it to win averages — leakage from "
            f"an unresolved peak fabricates residual power where the signal is "
            f"strongest."
        )
    elif verdict == "marginal":
        advice = (
            f"Be/Br = {ratio:.2f} at nperseg={nperseg}; the classical rule asks "
            f"for 0.25 or less. Raise nperseg to about {nperseg_res}, which "
            f"leaves {_welch_segment_count(n, nperseg_res, None)} segment(s) "
            f"here."
        )
    elif verdict == "too_few":
        fix, count = _suggest_nperseg(n, target_segments, nperseg, noverlap)
        advice = (
            f"{n_seg} segment(s) puts the coherence bias floor at "
            f"{floor:.2f}. nperseg={fix} would give {count}, but check the "
            f"resolution before shortening: "
            + (f"Be/Br would rise from {ratio:.2f} to "
               f"{resolution_bandwidth(fs, fix, window) / peak_bw:.2f}."
               if peak_bw else "no peak was found to check against.")
        )
    else:  # squeezed
        advice = (
            f"This record cannot satisfy both. Resolving the peak at "
            f"{peak_f:.4g} Hz needs nperseg>={nperseg_res}, which leaves "
            f"{_welch_segment_count(n, nperseg_res, None)} segment(s) of the "
            f"{n} samples available. About {duration_both:.4g} s of record "
            f"would give both, against the {n / fs:.4g} s you have. Until then "
            f"the choice is which error to prefer, and it should be made "
            f"knowingly."
        )

    return {
        "nperseg": nperseg,
        "n_segments": n_seg,
        "bias_floor": float(floor),
        "resolution_bw": be,
        "peak_freq": peak_f,
        "peak_bandwidth": peak_bw,
        "ratio": ratio,
        "resolved": resolved,
        "nperseg_for_resolution": nperseg_res,
        "duration_for_both": duration_both,
        "record_duration": float(n / fs),
        "verdict": verdict,
        "advice": advice,
    }

---

dspkit.spectral.autocorrelation(x, fs=None, normalize=True, max_lag=None)

Biased autocorrelation function (ACF) via FFT.

Uses the biased estimator (divides by N, not N-k) for better variance behaviour at large lags.

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	Input signal (zero-mean recommended; detrend first if needed).	required
fs	float or None	Sampling frequency [Hz]. If provided, the lag axis is in seconds; otherwise it is in samples.	None
normalize	bool	If True (default), normalise so that ACF[0] = 1.	True
max_lag	float or None	Maximum lag to return. Interpreted in seconds if fs is given, otherwise in samples. Defaults to the full one-sided ACF.	None

Returns:

Name	Type	Description
lags	ndarray	Lag axis (seconds if fs given, else samples).
acf	ndarray	Autocorrelation values.

Source code in dspkit/spectral.py

def autocorrelation(
    x: np.ndarray,
    fs: float | None = None,
    normalize: bool = True,
    max_lag: float | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Biased autocorrelation function (ACF) via FFT.

    Uses the biased estimator (divides by N, not N-k) for better
    variance behaviour at large lags.

    Parameters
    ----------
    x : array_like, shape (N,)
        Input signal (zero-mean recommended; detrend first if needed).
    fs : float or None
        Sampling frequency [Hz]. If provided, the lag axis is in seconds;
        otherwise it is in samples.
    normalize : bool
        If ``True`` (default), normalise so that ACF[0] = 1.
    max_lag : float or None
        Maximum lag to return. Interpreted in seconds if ``fs`` is given,
        otherwise in samples. Defaults to the full one-sided ACF.

    Returns
    -------
    lags : ndarray
        Lag axis (seconds if ``fs`` given, else samples).
    acf : ndarray
        Autocorrelation values.
    """
    x = np.asarray(x, dtype=float)
    N = len(x)

    # FFT-based circular correlation then truncate to causal (positive) lags
    Xf = np.fft.rfft(x, n=2 * N)
    acf_full = np.fft.irfft(Xf * np.conj(Xf))[:N]

    if normalize:
        acf_full = acf_full / acf_full[0]
    else:
        acf_full = acf_full / N

    lags = np.arange(N) / fs if fs is not None else np.arange(N, dtype=float)

    if max_lag is not None:
        cutoff = int(max_lag * fs) + 1 if fs is not None else int(max_lag) + 1
        lags = lags[:cutoff]
        acf_full = acf_full[:cutoff]

    return lags, acf_full

---

dspkit.spectral.cross_correlation(x, y, fs=None, normalize=True, max_lag=None)

Biased cross-correlation function (CCF) via FFT.

Computes the full two-sided CCF with lags from -(N-1) to +(N-1):

CCF[k] = (1/N) Σ_n x[n] · y[n + k]

A positive peak at lag k > 0 means y is the delayed copy: the pairing is x[n] with y[n+k], so if y[n] = x[n-d] the peak sits at k = +d and x leads y by d samples. (This docstring said the opposite until 2026-09-02; the formula above was always right. test_cross_correlation_lag_sign pins the direction.)

Parameters:

Name	Type	Description	Default
x	(array_like, shape(N))	Input signals. Must have the same length.	required
y	(array_like, shape(N))	Input signals. Must have the same length.	required
fs	float or None	Sampling frequency [Hz]. If provided, the lag axis is in seconds; otherwise it is in samples.	None
normalize	bool	If True (default), normalise so that max |CCF| ≤ 1. Specifically divides by sqrt(Rxx[0] · Ryy[0]), giving the cross-correlation coefficient — the same convention as numpy.corrcoef.	True
max_lag	float or None	Maximum absolute lag to return. Interpreted in seconds if fs is given, otherwise in samples.	None

Returns:

Name	Type	Description
lags	ndarray	Symmetric lag axis, running from -max_lag to +max_lag (seconds if fs given, else samples).
ccf	ndarray	Cross-correlation values.

Source code in dspkit/spectral.py

def cross_correlation(
    x: np.ndarray,
    y: np.ndarray,
    fs: float | None = None,
    normalize: bool = True,
    max_lag: float | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Biased cross-correlation function (CCF) via FFT.

    Computes the full two-sided CCF with lags from -(N-1) to +(N-1):

        CCF[k] = (1/N) Σ_n x[n] · y[n + k]

    A positive peak at lag k > 0 means **y is the delayed copy**: the pairing
    is x[n] with y[n+k], so if y[n] = x[n-d] the peak sits at k = +d and x
    leads y by d samples. (This docstring said the opposite until 2026-09-02;
    the formula above was always right. ``test_cross_correlation_lag_sign``
    pins the direction.)

    Parameters
    ----------
    x, y : array_like, shape (N,)
        Input signals. Must have the same length.
    fs : float or None
        Sampling frequency [Hz]. If provided, the lag axis is in seconds;
        otherwise it is in samples.
    normalize : bool
        If ``True`` (default), normalise so that ``max |CCF| ≤ 1``.
        Specifically divides by ``sqrt(Rxx[0] · Ryy[0])``, giving the
        cross-correlation coefficient — the same convention as
        ``numpy.corrcoef``.
    max_lag : float or None
        Maximum absolute lag to return. Interpreted in seconds if ``fs``
        is given, otherwise in samples.

    Returns
    -------
    lags : ndarray
        Symmetric lag axis, running from ``-max_lag`` to ``+max_lag``
        (seconds if ``fs`` given, else samples).
    ccf : ndarray
        Cross-correlation values.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    N = len(x)

    Xf = np.fft.fft(x, n=2 * N)
    Yf = np.fft.fft(y, n=2 * N)
    ccf_raw = np.fft.ifft(np.conj(Xf) * Yf).real

    # Rearrange to two-sided: lags -(N-1) … 0 … +(N-1)
    ccf_two = np.concatenate([ccf_raw[-(N - 1):], ccf_raw[:N]])

    if normalize:
        norm = np.sqrt(np.dot(x, x) * np.dot(y, y))
        if norm > 0:
            ccf_two = ccf_two / norm
    else:
        ccf_two = ccf_two / N

    lags_samples = np.arange(-(N - 1), N, dtype=float)
    lags = lags_samples / fs if fs is not None else lags_samples

    if max_lag is not None:
        cutoff = int(max_lag * fs) if fs is not None else int(max_lag)
        mask = np.abs(lags_samples) <= cutoff
        lags = lags[mask]
        ccf_two = ccf_two[mask]

    return lags, ccf_two