Learn more about System Identification Toolbox

spa - Estimate frequency response with fixed frequency resolution using spectral analysis

Syntax

G = spa(data) G = spa(data,winSize,freq) G = spa(data,winSize,freq,MaxSize) [G,phi,spectrum] = spa(data)

Description

G = spa(data) estimates frequency response (with uncertainty) and noise spectrum from time- or frequency-domain data. data is an iddata or idfrd object and can be complex valued. G is as an idfrd object. For time-series data, G is the estimated spectrum and standard deviation.

G = spa(data,winSize,freq) estimates frequency response at frequencies freq. freq is a row vector of values in rad/sec. winSize is a scalar integer that sets the size of the Hann window.

G = spa(data,winSize,freq,MaxSize) can improve computational performance using MaxSize to split the input-output data such that each segment contains fewer than MaxSize elements. MaxSize is a positive integer.

[G,phi,spectrum] = spa(data) estimates phi, which is the output disturbance and its uncertainty, and spetrum, which is the spectrum matrix for both the output and the input channels. For example, if z = [data.OutputData, data.InputData], spe contains the power spectrum estimate of z. G, phi, and spectrum are idfrd objects.

Definitions

Frequency Response Function

Frequency response function describes the steady-state response of a system to sinusoidal inputs. For a linear system, a sinusoidal input of a specific frequency results in an output that is also a sinusoid with the same frequency, but with a different amplitude and phase. The frequency response function describes the amplitude change and phase shift as a function of frequency.

To better understand the frequency response function, consider the following description of a linear, dynamic system:

where u(t) and y(t) are the input and output signals, respectively. G(q) is called the transfer function of the system—it captures the system dynamics that take the input to the output. The notation G(q)u(t) represents the following operation:

q is the shift operator, defined by the following equation:

G(q) is the frequency-response function, which is evaluated on the unit circle, G(q=e^iw).

Together, G(q=e^iw) and the output noise spectrum are the frequency-domain description of the system.

The frequency-response function estimated using the Blackman-Tukey approach is given by the following equation:

In this case, ^ represents approximate quantities. For a derivation of this equation, see the chapter on nonparametric time- and frequency-domain methods in System Identification: Theory for the User, Second Edition, by Lennart Ljung, Prentice Hall PTR, 1999.

Output Noise Spectrum

The output noise spectrum (spectrum of v(t)) is given by the following equation:

This equation for the noise spectrum is derived by assuming the linear relationship , that u(t) is independent of v(t), and the following relationships between the spectra:

where the noise spectrum is given by the following equation:

is the output-input cross-spectrum and is the input spectrum.

Alternatively, the disturbance v(t) can be described as filtered white noise:

where e(t) is the white noise with variance and the noise power spectrum is given by the following equation:

Examples

Estimate frequency response with fixed resolution at 128 equally spaced, logarithmic frequency values between 0 (excluded) and π:

g = spa(z); % z is an iddata object with Ts=1
bode(g)

Estimate frequency response with fixed resolution at logarithmically spaced frequencies:

% Define frequency vector
w = logspace(-2,pi,128);
% Compute frequency response
g= spa(z,[],w); % [] specifies the default lag window size
bode(g,'sd',3) % Bode plot of the transfer function
bode(g('noise'),'sd',3) % Noise spectrum
% Bode plots include confidence interval
% of 3 standard deviations

Algorithm

spa applies the Blackman-Tukey spectral analysis method by following these steps:

Computes the covariances and cross-covariance from u(t) and y(t):

This portion of the algorithm uses the covf function.
Computes the Fourier transforms of the covariances and the cross-covariance:

where is the Hann window with a width (lag size) of M. You can specify M to control the frequency resolution of the estimate, which is approximately equal 2π/M rad/sampling interval.
By default, this operation uses 128 equally spaced frequency values between 0 (excluded) and π, where w = [1:128]/128*pi/Ts and Ts is the sampling interval of that data set. The default lag size of the Hann window is M = min(length(data)/10,30). For default frequencies, uses fast Fourier transforms (FFT)—which is more efficient than for user-defined frequencies.
Compute the frequency-response function and the output noise spectrum .

spectrum is the spectrum matrix for both the output and the input channels. That is, if z = [data.OutputData, data.InputData], spectrum contains as spectrum data the matrix-valued power spectrum of z.