Documentation |
Time-domain correlation analysis refers to non-parametric estimation of the impulse response of dynamic systems as a finite impulse response (FIR) model from the data. Correlation analysis assumes a linear system and does not require a specific model structure.
Impulse response is the output signal that results when the input is an impulse and has the following definition for a discrete model:
$$\begin{array}{l}u(t)=0\text{}t0\\ u(t)=1\text{}t=0\end{array}$$
The response to an input u(t) is equal to the convolution of the impulse response, as follows:
$$y(t)={\displaystyle {\int}_{0}^{t}h\left(t-z\right)}\cdot u(z)dz$$
You can estimate impulse-response models from data with the following characteristics:
Real or complex data.
Single- or multiple-output data.
Time- or frequency-domain data with nonzero sampling time.
Time-domain data must be regularly sampled. You cannot use time-series data for correlation analysis.
Before you can perform this task, you must have:
Imported data into the System Identification app. See Import Time-Domain Data into the App. For supported data formats, see Data Supported by Correlation Analysis.
Performed any required data preprocessing operations. To improve the accuracy of your model, you should detrend your data. See Ways to Prepare Data for System Identification.
To estimate in the System Identification app using time-domain correlation analysis:
In the System Identification app, select Estimate > Correlation models to open the Correlation Model dialog box.
In the Time span (s) field, specify a scalar value as the time interval over which the impulse or step response is calculated. For a scalar time span T, the resulting response is plotted from -T/4 to T.
In the Order of whitening filter field, specify the filter order.
The prewhitening filter is determined by modeling the input as an autoregressive process of order N. The algorithm applies a filter of the form A(q)u(t)=u_F(t). That is, the input u(t) is subjected to an FIR filter A to produce the filtered signal u_F(t). Prewhitening the input by applying a whitening filter before estimation might improve the quality of the estimated impulse response g.
The order of the prewhitening filter, N, is the order of the A filter. N equals the number of lags. The default value of N is 10, which you can also specify as [].
In the Model Name field, enter the name of the correlation analysis model. The name of the model should be unique in the Model Board.
Click Estimate to add this model to the Model Board in the System Identification app.
In the Correlation Model dialog box, click Close.
Export the model to the MATLAB^{®} workspace for further analysis by dragging it to the To Workspace rectangle in the System Identification app.
View the transient response plot by selecting the Transient resp check box in the System Identification app. For more information about working with this plot and selecting to view impulse- versus step-response, see Impulse and Step Response Plots.
Before you can perform this task, you must have:
Input/output or frequency-response data. See Representing Time- and Frequency-Domain Data Using iddata Objects. For supported data formats, see Data Supported by Correlation Analysis.
Performed any required data preprocessing operations. If you use time-domain data, you can detrend it before estimation. See Ways to Prepare Data for System Identification.
Use impulseest to compute impulse response models. impulseest estimates a high-order, noncausal FIR model using correlation analysis. The resulting models are stored as idtf model objects and contain impulse-response coefficients in the model numerator.
To estimate the model m and plot the impulse or step response, use the following syntax:
m=impulseest(data,N); impulse(m,Time); step(m,Time);
where data is a single- or multiple-output iddata or idfrd object. N is a scalar value specifying the order of the FIR system corresponding to the time range 0:Ts:(N-1)*Ts, where Ts is the data sampling time.
You can also specify estimation options, such as regularizing kernel, pre-whitening filter order and data offsets, using impulseestOptions and pass them as an input to impulseest. For example:
opt = impulseestOptions('RegulKernel','TC')); m = impulseest(data,N,opt);
To view the confidence region for the estimated response, use impulseplot and stepplot to create the plot. Then use showConfidence.
For example:
h = stepplot(m,Time); showConfidence(h,3) % 3 std confidence region
Note: cra is an alternative method for computing impulse response from time-domain data only. |
Perform model analysis. See Validating Models After Estimation.
You can use impulse and step commands with output arguments to get the numerical impulse- and step-response vectors as a function of time, respectively.
To get the numerical response values:
Compute the FIR model by using impulseest, as described in Estimate Impulse-Response Models at the Command Line.
Apply the following syntax on the resulting model:
% To compute impulse-response data [y,t,~,ysd] = impulse(model) % To compute step-response data [y,t,~,ysd] = step(model)
where y is the response data, t is the time vector, and ysd is the standard deviations of the response.
You can use transient-response plots to estimate the input delay, or dead time, of linear systems. Input delay represents the time it takes for the output to respond to the input.
In the System Identification app. To view the transient response plot, select the Transient resp check box in the System Identification app. For example, the following step response plot shows a time delay of about 0.25 s before the system responds to the input.
Step Response Plot
At the command line. You can use impulseplot to plot the impulse response. The time delay is equal to the first positive peak in the transient response magnitude that is greater than the confidence region for positive time values.
For example, the following commands create an impulse-response plot with a 1-standard-deviation confidence region:
load dry2 ze = dry2(1:500); opt = impulseestOptions('RegulKernel','TC'); sys = impulseest(ze,40,opt); h = impulseplot(sys); showConfidence(h,1);
The resulting figure shows that the first positive peak of the response magnitude, which is greater than the confidence region for positive time values, occurs at 0.24 s.
Instead of using showConfidence, you can plot the confidence interval interactively, by right-clicking on the plot and selecting Characteristics > Confidence Region.
Correlation analysis refers to methods that estimate the impulse response of a linear model, without specific assumptions about model orders.
The impulse response, g, is the system's output when the input is an impulse signal. The output response to a general input, u(t), is obtained as the convolution with the impulse response. In continuous time:
$$y(t)={\displaystyle {\int}_{-\infty}^{t}g\left(\tau \right)u\left(t-\tau \right)}d\tau $$
In discrete-time:
$$y\left(t\right)={\displaystyle \sum _{k=1}^{\infty}g\left(k\right)u\left(t-k\right)}$$
The values of g(k) are the discrete time impulse response coefficients.
You can estimate the values from observed input-output data in several different ways. impulseest estimates the first n coefficients using the least-squares method to obtain a finite impulse response (FIR) model of order n.
Several important options are associated with the estimate:
Prewhitening — The input can be pre-whitened by applying an input-whitening filter of order PW to the data. This minimizes the effect of the neglected tail (k > n) of the impulse response.
A filter of order PW is applied such that it whitens the input signal u:
1/A = A(u)e, where A is a polynomial and e is white noise.
The inputs and outputs are filtered using the filter:
uf = Au, yf = Ay
The filtered signals uf and yf are used for estimation.
You can specify prewhitening using the PW name-value pair argument of impulseestOptions.
Regularization — The least-squares estimate can be regularized. This means that a prior estimate of the decay and mutual correlation among g(k) is formed and used to merge with the information about g from the observed data. This gives an estimate with less variance, at the price of some bias. You can choose one of the several kernels to encode the prior estimate.
This option is essential because, often, the model order n can be quite large. In cases where there is no regularization, n can be automatically decreased to secure a reasonable variance.
You can specify the regularizing kernel using the RegulKernel Name-Value pair argument of impulseestOptions.
Autoregressive Parameters — The basic underlying FIR model can be complemented by NA autoregressive parameters, making it an ARX model.
$$y\left(t\right)={\displaystyle \sum _{k=1}^{n}g\left(k\right)u\left(t-k\right)}-{\displaystyle \sum _{k=1}^{NA}{a}_{k}y\left(t-k\right)}$$
This gives both better results for small n and allows unbiased estimates when data are generated in closed loop. impulseest uses NA = 5 for t>0 and NA = 0 (no autoregressive component) for t<0.
Noncausal effects — Response for negative lags. It may happen that the data has been generated partly by output feedback:
$$u(t)={\displaystyle \sum _{k=0}^{\infty}h(k)y\left(t-k\right)}+r\left(t\right)$$
where h(k) is the impulse response of the regulator and r is a setpoint or disturbance term. The existence and character of such feedback h can be estimated in the same way as g, simply by trading places between y and u in the estimation call. Using impulseest with an indication of negative delays, $$\text{mi}=\text{impulseest}(data,nk,nb),\text{}nk0$$, returns a model mi with an impulse response
$$\left[h(-nk),h(-nk-1),\mathrm{...},h(0),g(1),g(2),\mathrm{...},g(nb+nk)\right]$$
aligned so that it corresponds to lags $$\left[nk,nk+1,\mathrm{..},0,1,2,\mathrm{...},nb+nk\right]$$. This is achieved because the input delay (InputDelay) of model mi is nk.
For a multi-input multi-output system, the impulse response g(k) is an ny-by-nu matrix, where ny is the number of outputs and nu is the number of inputs. The i–j element of the matrix g(k) describes the behavior of the ith output after an impulse in the jth input.