sidFreqBT

Python equivalent: sid.freq_bt

Estimate frequency response via Blackman-Tukey spectral analysis.

result = sidFreqBT(y, u)
result = sidFreqBT(y, [])
result = sidFreqBT(y, u, 'WindowSize', M, 'Frequencies', w)
result = sidFreqBT(y, u, M)
result = sidFreqBT(y, u, M, w)

Estimates the frequency response \(G(e^{jw})\) and noise spectrum from time-domain input/output data using the Blackman-Tukey method.

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

Inputs

Name	Description
y	Output data, (N x n_y) matrix. Column vector for SISO. For multiple trajectories: (N x n_y x L) array or cell array {y1, y2, ...} for variable-length data. Spectral estimates are ensemble-averaged across trajectories.
u	Input data, (N x n_u) matrix. Column vector for SISO. For multiple trajectories: (N x n_u x L) or cell array. Use [] for time series (output spectrum only).

Name-value options

Name	Description
'WindowSize'	Hann window lag size M. Default: min(floor(N/10), 30).
'Frequencies'	Frequency vector in rad/sample, in (0, pi]. Default: 128 linearly spaced values.
'SampleTime'	Sample time in seconds. Default: 1.0.

Outputs

Name	Description
result	Struct with fields:
.Frequency	(n_f x 1) frequency vector, rad/sample
.FrequencyHz	(n_f x 1) frequency vector, Hz
.Response	(n_f x n_y x n_u) complex frequency response
.ResponseStd	(n_f x n_y x n_u) standard deviation of Response
.NoiseSpectrum	(n_f x n_y x n_y) noise spectrum
.NoiseSpectrumStd	(n_f x n_y x n_y) standard deviation
.Coherence	(n_f x 1) squared coherence (SISO only)
.SampleTime	sample time in seconds
.WindowSize	window size M used
.DataLength	number of samples N
.NumTrajectories	number of trajectories L
.Method	'sidFreqBT'

Examples

% SISO system identification
N = 1000; u = randn(N,1);
y = filter([1], [1 -0.9], u) + 0.1*randn(N,1);
result = sidFreqBT(y, u);
sidBodePlot(result);

% Time series spectrum estimation
y = randn(500,1);
result = sidFreqBT(y, []);
sidSpectrumPlot(result);

% Custom window size and frequencies
w = linspace(0.01, pi, 256)';
result = sidFreqBT(y, u, 'WindowSize', 50, 'Frequencies', w);

% Multi-trajectory: average 5 independent experiments
L = 5; N = 1000;
y3d = zeros(N, 1, L); u3d = zeros(N, 1, L);
for l = 1:L
    u3d(:,1,l) = randn(N, 1);
    y3d(:,1,l) = filter([1], [1 -0.9], u3d(:,1,l)) + 0.1*randn(N,1);
end
result = sidFreqBT(y3d, u3d);

Algorithm

1. Compute biased sample covariances R_y, R_u, R_yu for lags 0..M
2. Apply Hann window and Fourier transform to obtain spectral estimates
3. Form G = \(Phi_yu\) / \(Phi_u\) and \(Phi_v\) = \(Phi_y\) - |\(Phi_yu\)|^2 / \(Phi_u\)
4. Compute asymptotic standard deviations

References

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

Specification

SPEC.md §2 — Blackman-Tukey Spectral Analysis

See also

sidFreqBTFDR, sidFreqETFE, sidBodePlot, sidSpectrumPlot

Changelog

• 2026-03-24: First version by Pedro Lourenço.