spectrum

Output power spectrum of time series models

collapse all in page

Syntax

spectrum(sys)
spectrum(sys,{wmin, wmax})
spectrum(sys,w)
spectrum(sys1,...,sysN,w)
ps = spectrum(sys,w)
[ps,w] = spectrum(sys)
[ps,w,sdps] = spectrum(sys)

Description

spectrum(sys) creates an output power spectrum plot of the identified time series model sys. The frequency range and number of points are chosen automatically.

sys is a time series model, which represents the system:

y(t)=He(t)

Where, e(t) is a Gaussian white noise and y(t) is the observed output.

spectrum plots abs(H'H), scaled by the variance of e(t) and the sample time.

If sys is an input-output model, it represents the system:

y(t)=Gu(t)+He(t)

Where, u(t) is the measured input, e(t) is a Gaussian white noise and y(t) is the observed output.

In this case, spectrum plots the spectrum of the disturbance component He(t).

spectrum(sys,{wmin, wmax}) creates a spectrum plot for frequencies ranging from wmin to wmax.

spectrum(sys,w) creates a spectrum plot using the frequencies specified in the vector w.

spectrum(sys1,...,sysN,w) creates a spectrum plot of several identified models on a single plot. The w argument is optional.

You can specify a color, line style and marker for each model. For example:

spectrum(sys1,'r',sys2,'y--',sys3,'gx');

ps = spectrum(sys,w) returns the power spectrum amplitude of sys for the specified frequencies, w. No plot is drawn on the screen.

[ps,w] = spectrum(sys) returns the frequency vector, w, for which the output power spectrum is plotted.

[ps,w,sdps] = spectrum(sys) returns the estimated standard deviations of the power spectrum.

For discrete-time models with sample time Ts, spectrum uses the transformation z = exp(j*w*Ts) to map the unit circle to the real frequency axis. The spectrum is only plotted for frequencies smaller than the Nyquist frequency pi/Ts, and the default value 1 (time unit) is assumed when Ts is unspecified.

Input Arguments

sys

Identified model.

If sys is a time series model, it represents the system:

y(t)=He(t)

Where, e(t) is a Gaussian white noise and y(t) is the observed output.

If sys is an input-output model, it represents the system:

y(t)=Gu(t)+He(t)

Where, u(t) is the measured input, e(t) is a Gaussian white noise and y(t) is the observed output.

wmin

Minimum frequency of the frequency range for which the output power spectrum is plotted.

Specify wmin in rad/TimeUnit, where TimeUnit is sys.TimeUnit.

wmax

Maximum frequency of the frequency range for which the output power spectrum is plotted.

Specify wmax in rad/TimeUnit, where TimeUnit is sys.TimeUnit.

w

Frequencies for which the output power spectrum is plotted.

Specify w in rad/TimeUnit, where TimeUnit is sys.TimeUnit.

sys1,...,sysN

Identified systems for which the output power spectrum is plotted.

Output Arguments

ps

Power spectrum amplitude.

If sys has Ny outputs, then ps is an array of size [Ny Ny length(w)]. Where ps(:,:,k) corresponds to the power spectrum for the frequency at w(k).

For amplitude values in dB, type psdb = 10*log10(ps).

w

Frequency vector for which the output power spectrum is plotted.

sdps

Estimated standard deviation of the power spectrum.

Examples

collapse all

Open Live Script

Load the estimation data. 

load iddata1 z1;

Estimate a single-input single-output state-space model. 

sys = n4sid(z1,2);

Plot the noise spectrum for the model.

spectrum(sys);

Open Live Script

Load the time-series estimation data. 

load iddata9 z9

Estimate a fourth-order AR model using a least-squares approach.

sys = ar(z9,4,'ls');

Plot the output spectrum of the model. 

spectrum(sys);

Open Live Script

Create an input consisting of five sinusoids spread over the whole frequency interval. Compare the spectrum of this signal with that of its square. The frequency splitting (the square having spectral support at other frequencies) reveals the nonlinearity involved.

u = idinput([100 1 20],'sine',[],[],[5 10 1]);
u = iddata([],u,1,'per',100);
u2 = u.u.^2;
u2 = iddata([],u2,1,'per',100);
spectrum(etfe(u),'r*',etfe(u2),'+')

See Also

| | | | | |

Introduced in R2012a

System Identification Toolbox Documentation

Support