freq_bt¶

MATLAB equivalent: sidFreqBT

freq_bt ¶

Blackman-Tukey spectral analysis for frequency response estimation.

freq_bt ¶

freq_bt(y: ndarray | list, u: ndarray | list | None = None, *, window_size: int | None = None, frequencies: ndarray | None = None, sample_time: float = 1.0) -> FreqResult

Estimate frequency response via Blackman-Tukey spectral analysis.

Estimates the frequency response G(e^{jw}) and noise power spectrum from time-domain input/output data using the Blackman-Tukey method. The algorithm computes biased sample covariances, applies a Hann lag window, and Fourier-transforms to obtain spectral estimates.

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

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`
`window_size`	`int`	Hann lag window size M. Must satisfy `M >= 2`. If M exceeds `N // 2` it is silently reduced. Default is `min(N // 10, 30)`.	`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 (int) -- Lag window size M used.
data_length (int) -- Number of samples N per trajectory.
num_trajectories (int) -- Number of trajectories.
method (str) -- 'freq_bt'.

Raises:

Type	Description
`SidError`	If sample_time is not positive (code: `'bad_ts'`).
`SidError`	If window_size is less than 2 (code: `'bad_window_size'`).
`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_bt(y, u)
>>> result.response.shape
(128,)

Time-series spectrum estimation:

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

Custom window size and frequencies:

>>> w = np.linspace(0.01, np.pi, 256)
>>> result = sid.freq_bt(y, u, window_size=50, 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_bt(y3d, u3d)

Notes

Algorithm:

Validate and orient input data; detect time-series vs. SISO vs. MIMO.
Compute biased sample covariances R_yy, R_uu, R_yu for lags 0 .. M (ensemble-averaged for multi-trajectory data).
Apply Hann lag window and Fourier-transform to obtain spectral estimates Phi_y, Phi_u, Phi_yu. Uses the FFT fast path when the frequency grid matches the default 128-point linear grid, and direct DFT otherwise.
Form the transfer function estimate G = Phi_yu / Phi_u (SISO) or G = Phi_yu * Phi_u^{-1} (MIMO), and the noise spectrum Phi_v = Phi_y - |Phi_yu|^2 / Phi_u (SISO) or Phi_v = Phi_y - Phi_yu * Phi_u^{-1} * Phi_yu' (MIMO).
Compute asymptotic standard deviations using Ljung (1999) formulae.

Specification: SPEC.md S2 -- Blackman-Tukey Spectral Analysis

References

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