spa (System Identification Toolbox)

spa	Examples See Also

Estimate frequency response and spectrum by spectral analysis.

Syntax

[g,phiv] = spa(z)
[g,phiv,z_spe] = spa(z,M,w,maxsize,T)

Description

spa estimates the transfer function g and the noise spectrum

phiv =

of the general linear model

where

is the spectrum of

(t).

Matrix z contains the output-input data z = [y u], where y and u are column vectors. If there are several inputs, u is a matrix with one column for each input. The data may be complex-valued.

g is returned in the frequency function format (see freqfunc) with the estimate of

at the frequencies

specified by row vector w. The default value of w is

w = [1:128]/128*pi/T

phiv is returned with the autospectrum estimate of

at the same frequencies. Both outputs are returned with estimated standard deviations (see freqfunc).

M is the length of the lag window used in the calculations. The default value is

M = min(30,length(z)/10)

Changing the value of M exchanges bias for variance in the spectral estimate. As M is increased, the estimated functions show more detail, but are more corrupted by noise. The sharper peaks a true frequency function has, the higher M it needs. See etfe as an alternative for narrowband signals and systems.

T is the sampling interval and maxsize controls the memory-speed trade-off (see auxvar).

For time series z = y, g is returned with the estimated output spectrum and its estimated standard deviation.

Important: For multi-output data the argument M must be entered as a row vector of the same length as the number of outputs. This is the way the distinction between inputs and outputs in z is clarified. For default window size use in the multi-output case

M = [-1, -1, ..., -1]

The optional third output argument z_spe gives directly the spectrum matrix of z as follows:

reshape(z_spe(:,k),nz,nz) = The spectrum S at frequency W(k)

where nz is the number of channels in the data matrix z and

Here win(m) is weight at lag m of an M-size Hamming window and W(k) is the frequency value i rad/s. Note that ' denotes complex-conjugate transpose.

The normalization of the spectrum differs from the one used by spectrum in the Signal Processing Toolbox . See "Some Special Topics" on page 3-68 in the User's Guide for a more precise definition.

Examples

With default frequencies

g = spa(z);
bodeplot(g)

With logarithmically spaced frequencies

w = logspace(-2,pi,128);
[g,phiv] = spa(z,[],w); 
% (empty matrix gives default)
bodeplot([g phiv],3)

plots the estimated spectrum together with confidence intervals of three standard deviations.

Algorithm

The covariance function estimates are computed using covf. These are multiplied by a Hamming window of lag size M and then Fourier transformed. The relevant ratios and differences are then formed. For the default frequencies, this is done using FFT, which is more efficient than for user-defined frequencies. For multi-variable systems, a straightforward for-loop is used.

Note that M =

is in Table 6.1 of Ljung (1987). The standard deviations are computed as on pages 155-156 and 264 of Ljung (1987).