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)

Description

G = spa(data) estimates the frequency response, along with uncertainty , and the noise spectrum from time- or frequency-domain data data. If data is a time series, spa(data) returns the output power spectrum along with uncertainty. spa computes the spectra at 128 equally spaced frequency values between 0 (excluded) and π, using a Hann window.

G = spa(data,winSize,freq) estimates the frequency response at the frequencies contained in freq, using a Hann window with size winSize.

G = spa(data,winSize,freq,maxSize) splits the input-output data into segments, each segment containing fewer than maxSize elements. Use this syntax to improve computational performance.

Examples

Estimate Frequency Response

Estimate frequency response with fixed resolution at 128 equally spaced, logarithmic frequency values between 0 (excluded) and $$\pi$$.

load iddata3; 
g = spa(z3); 
bode(g)

Estimate Frequency Response at Specified Frequencies

Define the frequency vector.

w = logspace(-2,pi,128);

Compute the frequency response.

load iddata3;
g = spa(z3,[],w);

[] specifies the default lag window size.

Plot the Bode response and disturbance spectrum with confidence interval of 3 standard deviations.

h = bodeplot(g);
showConfidence(h,3)

figure
h = spectrumplot(g);
showConfidence(h,3)

Input Arguments

data — Input-output data

iddata object | idfrd object

Input-output data, specified as an iddata object or an idfrd object that can contain complex values. data can also contain time series data with only output.

winSize — Window size

scalar integer

Hann window size, also known as lag size, specified as a scalar integer. By default, the function sets the window size to min(length(data)/10,30).

freq — Frequencies

row vector

Frequencies at which to estimate spectral response, specified as a row vector in units of rad/timeUnit, where timeUnit refers to the timeUnit property of data. By default, the function sets freq to a vector of 128 values in the range (0,π], evenly spaced logarithmically. For discrete-time models, set freq within the Nyquist frequency bound.

maxSize — Maximum segment size

250e3 | positive integer

Maximum size of segments within data, specified as a positive integer. If you omit this argument, the function performs estimation using the full data set in data rather than segmenting the data.

Output Arguments

G — Frequency response and noise spectrum

idfrd object

Frequency response with uncertainty and noise spectrum, specified as an idfrd object. For time series data, G is the estimated spectrum and standard deviation.

Information about the estimation results and options used is stored in the Report property of the model. Report has the following fields.

Report Field	Description
Status	Summary of the model status, which indicates whether the model was created by construction or obtained by estimation.
Method	Estimation command used.
windowSize	Size of the Hann window.
DataUsed	Attributes of the data used for estimation. Structure with the following fields: Field Description Name Name of the data set. Type Data type. Length Number of data samples. Ts Sample time. This is equivalent to Data.Ts. InterSample Input intersample behavior. One of the following values: 'zoh' — Zero-order hold maintains a piecewise-constant input signal between samples. 'foh' — First-order hold maintains a piecewise-linear input signal between samples. 'bl' — Band-limited behavior specifies that the continuous-time input signal has zero power above the Nyquist frequency. The value of Intersample has no effect on estimation results for discrete-time models. InputOffset Offset removed from time-domain input data during estimation. OutputOffset Offset removed from time-domain output data during estimation.

More About

Frequency Response Function

A 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:

$$y(t)=G(q)u(t)+v(t)$$

Here 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:

$$G(q)u(t)={\displaystyle \sum _{k=1}^{\infty }g(k)u(t-k)}$$

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

$$G(q)={\displaystyle \sum _{k=1}^{\infty }g(k)q^{-k}}{q}^{-1}u(t)=u(t-1)$$

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

Together, G(q=e^iω) and the output noise spectrum $${\widehat{\Phi }}_v(\omega )$$ are the frequency-domain description of the system.

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

$${\widehat{G}}_N\left(e^{i\omega }\right)=\frac{{\widehat{\Phi }}_{yu}\left(\omega \right)}{{\widehat{\Phi }}_u\left(\omega \right)}$$

In this case, ^ represents approximate quantities. For a derivation of this equation, see the chapter on nonparametric time- and frequency-domain methods in [1].

Output Noise Spectrum

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

$${\widehat{\Phi }}_v\left(\omega \right)={\widehat{\Phi }}_y\left(\omega \right)-\frac{{\left|{\widehat{\Phi }}_{yu}\left(\omega \right)\right|}^2}{{\widehat{\Phi }}_u\left(\omega \right)}$$

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

$${\Phi }_y(\omega )={\left|G\left(e^{i\omega }\right)\right|}^2{\Phi }_u(\omega )+{\Phi }_v(\omega )$$

$${\Phi }_{yu}(\omega )=G\left(e^{i\omega }\right){\Phi }_u(\omega )$$

Here, the noise spectrum is given by the following equation:

$${\Phi }_v(\omega )\equiv {\displaystyle \sum _{\tau =-\infty }^{\infty }{R}_v}(\tau )e^{-iw\tau }$$

$${\widehat{\Phi }}_{yu}(\omega )$$ is the output-input cross-spectrum and $${\widehat{\Phi }}_u(\omega )$$ is the input spectrum.

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

$$v(t)=H(q)e(t)$$

Here, e(t) is the white noise with variance $$\lambda$$ and the noise power spectrum is given by the following equation:

$${\Phi }_v(\omega )=\lambda {\left|H\left(e^{i\omega }\right)\right|}^2$$

Algorithms

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

1. Compute the covariances and cross-covariance from u(t) and y(t):

   $$\begin{array}{l}{\widehat{R}}_y\left(\tau \right)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N}y\left(t+\tau \right)y\left(t\right)}\\ {\widehat{R}}_u\left(\tau \right)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N}u\left(t+\tau \right)u\left(t\right)}\\ {\widehat{R}}_{yu}\left(\tau \right)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N}y\left(t+\tau \right)u\left(t\right)}\end{array}$$

2. Compute the Fourier transforms of the covariances and the cross-covariance:

   $$\begin{array}{l}{\widehat{\Phi }}_y(\omega )={\displaystyle \sum _{\tau =-M}^M{\widehat{R}}_y}(\tau )W_M(\tau )e^{-i\omega \tau }\\ {\widehat{\Phi }}_u(\omega )={\displaystyle \sum _{\tau =-M}^M{\widehat{R}}_u}(\tau )W_M(\tau )e^{-i\omega \tau }\\ {\widehat{\Phi }}_{yu}(\omega )={\displaystyle \sum _{\tau =-M}^M{\widehat{R}}_{yu}}(\tau )W_M(\tau )e^{-i\omega \tau }\end{array}$$

   where $$W_M(\tau )$$ 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/sample time.

   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 sample time of that data set. The default lag size of the Hann window is M = min(length(data)/10,30). For default frequencies, the operation uses fast Fourier transforms (FFT), which is more efficient than for user-defined frequencies.

   Note

   M =γ is in Table 6.1 of [1]. Standard deviations are on pages 184 and 188 in [1].

3. Compute the frequency-response function $${\widehat{G}}_N\left(e^{i\omega }\right)$$ and the output noise spectrum $${\widehat{\Phi }}_v(\omega )$$.

   $${\widehat{G}}_N\left(e^{i\omega }\right)=\frac{{\widehat{\Phi }}_{yu}\left(\omega \right)}{{\widehat{\Phi }}_u\left(\omega \right)}$$

   $${\Phi }_v(\omega )\equiv {\displaystyle \sum _{\tau =-\infty }^{\infty }{R}_v}(\tau )e^{-iw\tau }$$

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.

$$S={\displaystyle \sum _{m=-M}^{M}Ez\left(t+m\right)z}{\left(t\right)}^{\prime }{W}_M\left(T_s\right)\exp \left(-i\omega m\right)$$

' is a complex-conjugate transpose.

References

[1] Ljung, Lennart. System Identification: Theory for the User, Second Ed., Prentice Hall PTR, 1999.