freq_btfdr¶

MATLAB equivalent: sidFreqBTFDR

freq_btfdr ¶

Blackman-Tukey spectral analysis with frequency-dependent resolution.

freq_btfdr ¶

freq_btfdr(y: ndarray | list, u: ndarray | list | None = None, *, resolution: float | ndarray | None = None, frequencies: ndarray | None = None, sample_time: float = 1.0) -> FreqResult

Estimate frequency response via Blackman-Tukey with frequency-dependent resolution.

Like :func:sid.freq_bt, but the Hann lag-window size varies across frequencies. The user specifies a resolution parameter (in rad/sample) instead of a fixed window size. Finer resolution (smaller value) uses a larger window and gives lower variance but coarser frequency detail, while coarser resolution (larger value) uses a smaller window.

This is an open-source replacement for the System Identification Toolbox function spafdr.

Parameters:

Name	Type	Description	Default
`y`	`ndarray or list of ndarray, shape (N, ny) or (N, ny, L)`	Output data. A 1-D array is treated as a single-channel signal. For multiple trajectories pass a 3-D array `(N, ny, L)` or a list of 1-D / 2-D arrays (variable-length trajectories are trimmed to the shortest). Spectral estimates are ensemble-averaged across trajectories.	required
`u`	`ndarray, list of ndarray, or None`	Input data. Same shape conventions as y. Pass `None` for time-series mode (output spectrum only). Default is `None`.	`None`
`resolution`	`float, ndarray, or None`	Frequency resolution in rad/sample. A scalar value applies uniformly to all frequencies; a vector of the same length as frequencies sets per-frequency resolution. Must be positive. Default is `2 * pi / min(N // 10, 30)` (clamped so the implied window size is at least 2).	`None`
`frequencies`	`(ndarray, shape(nf))`	Frequency vector in rad/sample. All values must lie in the interval `(0, pi]`. Default is 128 linearly spaced values `k * pi / 128` for `k = 1, ..., 128`.	`None`
`sample_time`	`float`	Sample time in seconds. Must be positive. Default is `1.0`.	`1.0`

Returns:

Type Description

FreqResult

Frozen dataclass with fields:

frequency (ndarray, shape (nf,)) -- Frequency vector, rad/sample.
frequency_hz (ndarray, shape (nf,)) -- Frequency vector, Hz.
response (ndarray or None) -- Complex frequency response. Shape (nf,) for SISO, (nf, ny, nu) for MIMO, or None in time-series mode.
response_std (ndarray or None) -- Standard deviation of response, same shape. None in time-series mode.
noise_spectrum (ndarray) -- Noise power spectrum (or output spectrum in time-series mode). Shape (nf,) for SISO / time-series, (nf, ny, ny) for MIMO.
noise_spectrum_std (ndarray) -- Standard deviation of noise_spectrum, same shape.
coherence (ndarray or None) -- Squared coherence, shape (nf,). SISO only; None for MIMO or time-series.
sample_time (float) -- Sample time in seconds.
window_size (ndarray, shape (nf,)) -- Per-frequency lag window sizes M_k.
data_length (int) -- Number of samples N per trajectory.
num_trajectories (int) -- Number of trajectories.
method (str) -- 'freq_btfdr'.

Raises:

Type	Description
`SidError`	If sample_time is not positive (code: `'bad_ts'`).
`SidError`	If any resolution value is not positive (code: `'bad_resolution'`).
`SidError`	If resolution is a vector whose length does not match the frequency grid (code: `'bad_resolution'`).
`SidError`	If any frequency is outside `(0, pi]` (code: `'bad_freqs'`).
`SidError`	If data contains NaN/Inf (code: `'non_finite'`), is complex (code: `'complex_data'`), or is too short (code: `'too_short'`). These are raised by the data validation layer.

Examples:

Basic SISO system identification:

>>> import numpy as np
>>> import sid
>>> N = 1000; rng = np.random.default_rng(0)
>>> u = rng.standard_normal(N)
>>> from scipy.signal import lfilter
>>> y = lfilter([1], [1, -0.9], u) + 0.1 * rng.standard_normal(N)
>>> result = sid.freq_btfdr(y, u, resolution=0.3)
>>> result.response.shape
(128,)

Time-series spectrum estimation:

>>> y = rng.standard_normal(500)
>>> result = sid.freq_btfdr(y)
>>> result.noise_spectrum.shape
(128,)

Per-frequency resolution vector:

>>> w = np.linspace(0.01, np.pi, 64)
>>> R = np.linspace(0.1, 1.0, 64)
>>> result = sid.freq_btfdr(y, u, resolution=R, frequencies=w)

Multi-trajectory (ensemble-averaged):

>>> L = 5; N = 1000
>>> u3d = rng.standard_normal((N, 1, L))
>>> y3d = np.zeros_like(u3d)
>>> for l in range(L):
...     y3d[:, 0, l] = lfilter([1], [1, -0.9], u3d[:, 0, l]) + 0.1 * rng.standard_normal(N)
>>> result = sid.freq_btfdr(y3d, u3d, resolution=0.3)

Notes

Algorithm:

Validate and orient input data; detect time-series vs. SISO vs. MIMO.
Convert the resolution parameter to per-frequency window sizes M_k = ceil(2 * pi / R_k), clamped to [2, N // 2].
Pre-compute biased sample covariances up to max(M_k) using ensemble averaging for multi-trajectory data.
For each frequency w_k:

a. Truncate covariances to lag M_k and compute a local Hann window of size M_k. b. Compute windowed spectral estimates via a direct single-frequency DFT (not FFT, since the window size varies per frequency). c. Form G, noise spectrum, coherence, and asymptotic uncertainty using the local window norm C_W.

The per-frequency approach gives smooth bias-variance trade-offs across the frequency axis, at the cost of a slower per-frequency loop compared to the fixed-window :func:sid.freq_bt.

Specification: SPEC.md S5 -- Frequency-Dependent Resolution

References

.. [1] Ljung, L., "System Identification: Theory for the User", 2nd ed., Prentice Hall, 1999. Sections 6.3--6.4.