Documentation |
You can estimate frequency-response models and visualize the responses on a Bode plot, which shows the amplitude change and the phase shift as a function of the sinusoid frequency.
The 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.
For a discrete-time system sampled with a time interval T, the frequency-response model G(z) relates the Z-transforms of the input U(z) and output Y(z):
$$Y(z)=G(z)U(z)$$
In other words, the frequency-response function, G(e^{iwT}), is the Laplace transform of the impulse response that is evaluated on the imaginary axis. The frequency-response function is the transfer function G(z) evaluated on the unit circle.
The estimation result is an idfrd model, which stores the estimated frequency response and its covariance.
You can estimate spectral analysis models from data with the following characteristics:
Complex or real data.
Time- or frequency-domain iddata or idfrd data object. To learn more about estimating time-series models, see Time-Series Model Identification.
Single- or multiple-output data.
You must have already imported your data into the app and performed any necessary preprocessing operations. For more information, see Data Preparation.
To estimate frequency-response models in the System Identification app:
In the System Identification app, select Estimate > Spectral models to open the Spectral Model dialog box.
In the Method list, select the spectral analysis method you want to use. For information about each method, see Selecting the Method for Computing Spectral Models.
Specify the frequencies at which to compute the spectral model in one of the following ways:
In the Frequencies field, enter either a vector of values, a MATLAB^{®} expression that evaluates to a vector, or a variable name of a vector in the MATLAB workspace. For example, logspace(-1,2,500).
Use the combination of Frequency Spacing and Frequencies to construct the frequency vector of values:
In the Frequency Spacing list, select Linear or Logarithmic frequency spacing.
In the Frequencies field, enter the number of frequency points.
For time-domain data, the frequency ranges from 0 to the Nyquist frequency. For frequency-domain data, the frequency ranges from the smallest to the largest frequency in the data set.
In the Frequency Resolution field, enter the frequency resolution, as described in Controlling Frequency Resolution of Spectral Models . To use the default value, enter default or, equivalently, the empty matrix [].
In the Model Name field, enter the name of the correlation analysis model. The model name should be unique in the Model Board.
Click Estimate to add this model to the Model Board in the System Identification app.
In the Spectral Model dialog box, click Close.
To view the frequency-response plot, select the Frequency resp check box in the System Identification app. For more information about working with this plot, see Frequency Response Plots.
To view the estimated disturbance spectrum, select the Noise spectrum check box in the System Identification app. For more information about working with this plot, see Noise Spectrum Plots.
Validate the model after estimating it. For more information, see Model Validation.
To export the model to the MATLAB workspace, drag it to the To Workspace rectangle in the System Identification app. You can retrieve the responses from the resulting idfrd model object using the bode or nyquist command.
You can use the etfe, spa, and spafdr commands to estimate spectral models. The following table provides a brief description of each command and usage examples.
The resulting models are stored as idfrd model objects. For detailed information about the commands and their arguments, see the corresponding reference page.
Commands for Frequency Response
Command | Description | Usage |
---|---|---|
etfe | Estimates an empirical transfer function using Fourier analysis. | To estimate a model m, use the following syntax: m=etfe(data) |
spa | Estimates a frequency response with a fixed frequency resolution using spectral analysis. | To estimate a model m, use the following syntax: m=spa(data) |
spafdr | Estimates a frequency response with a variable frequency resolution using spectral analysis. | To estimate a model m, use the following syntax: m=spafdr(data,R,w) where R is the resolution vector and w is the frequency vector. |
Validate the model after estimating it. For more information, see Model Validation.
This section describes how to select the method for computing spectral models in the estimation procedures Estimate Frequency-Response Models in the App and Estimate Frequency-Response Models at the Command Line.
You can choose from the following three spectral-analysis methods:
etfe (Empirical Transfer Function Estimate)
For input-output data. This method computes the ratio of the Fourier transform of the output to the Fourier transform of the input.
For time-series data. This method computes a periodogram as the normalized absolute squares of the Fourier transform of the time series.
ETFE works well for highly resonant systems or narrowband systems. The drawback of this method is that it requires linearly spaced frequency values, does not estimate the disturbance spectrum, and does not provide confidence intervals. ETFE also works well for periodic inputs and computes exact estimates at multiples of the fundamental frequency of the input and their ratio.
spa (SPectral Analysis)
This method is the Blackman-Tukey spectral analysis method, where windowed versions of the covariance functions are Fourier transformed.
spafdr (SPectral Analysis with Frequency Dependent Resolution)
This method is a variant of the Blackman-Tukey spectral analysis method with frequency-dependent resolution. First, the algorithm computes Fourier transforms of the inputs and outputs. Next, the products of the transformed inputs and outputs with the conjugate input transform are smoothed over local frequency regions. The widths of the local frequency regions can vary as a function of frequency. The ratio of these averages computes the frequency-response estimate.
This section supports the estimation procedures Estimate Frequency-Response Models in the App and Estimate Frequency-Response Models at the Command Line.
Frequency resolution is the size of the smallest frequency for which details in the frequency response and the spectrum can be resolved by the estimate. A resolution of 0.1 rad/s means that the frequency response variations at frequency intervals at or below 0.1 rad/s are not resolved.
Specifying the frequency resolution for etfe and spa is different than for spafdr.
For etfe and spa, the frequency resolution is approximately equal to the following value:
$$\frac{2\pi}{M}\left(\frac{\text{radians}}{\text{samplinginterval}}\right)$$
M is a scalar integer that sets the size of the lag window. The value of M controls the trade-off between bias and variance in the spectral estimate.
The default value of M for spa is good for systems without sharp resonances. For etfe, the default value of M gives the maximum resolution.
A large value of M gives good resolution, but results in more uncertain estimates. If a true frequency function has sharp peak, you should specify higher M values.
In case of etfe and spa, the frequency response is defined over a uniform frequency range, 0-F_{s}/2 radians per second, where F_{s} is the sampling frequency—equal to twice the Nyquist frequency. In contrast, spafdr lets you increase the resolution in a specific frequency range, such as near a resonance frequency. Conversely, you can make the frequency grid coarser in the region where the noise dominates—at higher frequencies, for example. Such customizing of the frequency grid assists in the estimation process by achieving high fidelity in the frequency range of interest.
For spafdr, the frequency resolution around the frequency k is the value R(k). You can enter R(k) in any one of the following ways:
Scalar value of the constant frequency resolution value in radians per second.
Vector of frequency values the same size as the frequency vector.
Expression using MATLAB workspace variables and evaluates to a resolution vector that is the same size as the frequency vector.
The default value of the resolution for spafdr is twice the difference between neighboring frequencies in the frequency vector.
If the input data is marked as periodic and contains an integer number of periods (data.Period is an integer), etfe computes the frequency response at frequencies $${\scriptscriptstyle \frac{2\pi k}{T}}\left({\scriptscriptstyle \frac{\text{k}}{\text{Period}}}\right)\text{where}k=1,2,\mathrm{...},\text{Period}$$.
For periodic data, the frequency resolution is ignored.
The spectrum of a signal is the square of the Fourier transform of the signal. The spectral estimate using the commands spa, spafdr, and etfe is normalized by the sampling interval T:
$${\Phi}_{y}(\omega )=T{\displaystyle \sum _{k=-M}^{M}{R}_{y}}(kT){e}^{-iwT}{W}_{M}(k)$$
where W_{M}(k) is the lag window, and M is the width of the lag window. The output covariance R_{y}(kT) is given by the following discrete representation:
$${\widehat{R}}_{y}(kT)=\frac{1}{N}{\displaystyle \sum _{l=1}^{N}y(lT-kT)y(lT)}$$
Because there is no scaling in a discrete Fourier transform of a vector, the purpose of T is to relate the discrete transform of a vector to the physically meaningful transform of the measured signal. This normalization sets the units of $${\Phi}_{y}(\omega )$$ as power per radians per unit time, and makes the frequency units radians per unit time.
The scaling factor of T is necessary to preserve the energy density of the spectrum after interpolation or decimation.
By Parseval's theorem, the average energy of the signal must equal the average energy in the estimated spectrum, as follows:
$$\begin{array}{l}E{y}^{2}(t)=\frac{1}{2\pi}{\displaystyle {\int}_{-\pi /T}^{\pi /T}{\Phi}_{y}}(\omega )d\omega \\ S1\equiv E{y}^{2}(t)\\ S2\equiv \frac{1}{2\pi}{\displaystyle {\int}_{-\pi /T}^{\pi /T}{\Phi}_{y}}(\omega )d\omega \end{array}$$
To compare the left side of the equation (S1) to the right side (S2), enter the following commands. In this code, phiy contains $${\Phi}_{y}(\omega )$$ between $$\omega =0$$ and $$\omega ={\scriptscriptstyle \raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$T$}\right.}$$ with the frequency step given as follows:
$$\left(\frac{\pi}{T\cdot \text{length(phiy)}}\right)$$
load iddata1 % Create a time-series iddata object. y = z1(:,1,[]); % Define sample interval from the data. T = y.Ts; % Estimate the frequency response. sp = spa(y); % Remove spurious dimensions phiy = squeeze(sp.spec); % Compute average energy from the estimated energy spectrum, where S1 is % scaled by T. S1 = sum(phiy)/length(phiy)/T % Compute average energy of the signal. S2 = sum(y.y.^2)/size(y,1)
S1 = 19.2076 S2 = 19.4646
Thus, the average energy of the signal approximately equals the average energy in the estimated spectrum.