Estimate parameters of ARX or AR models recursively
thm = rarx(z,nn,adm,adg)
[thm,yhat,P,phi] = rarx(z,nn,adm,adg,th0,P0,phi0)
thm = rarx(z,nn,adm,adg) estimates the parameters thm of single-output ARX model from input-output data z and model orders nn using the algorithm specified by adm and adg. If z is a time series y and nn = na, rarx estimates the parameters of a single-output AR model.
[thm,yhat,P,phi] = rarx(z,nn,adm,adg,th0,P0,phi0) estimates the parameters thm, the predicted output yhat, final values of the scaled covariance matrix of the parameters P, and final values of the data vector phi of single-output ARX model from input-output data z and model orders nn using the algorithm specified by adm and adg. If z is a time series y and nn = na, rarx estimates the parameters of a single-output AR model.
The general ARX model structure is:
The orders of the ARX model are:
Models with several inputs are defined, as follows:
A(q)y(t) = B1(q)u1(t–nk1)+...+Bnuunu(t–nknu)+e(t)
Name of the matrix iddata object that represents the input-output data or a matrix z = [y u], where y and u are column vectors.
For multiple-input models, the u matrix contains each input as a column vector:
u = [u1 ... unu]
For input-output models, specifies the structure of the ARX model as:
nn = [na nb nk]
where na and nb are the orders of the ARX model, and nk is the delay.
For multiple-input models, nb and nk are row vectors that define orders and delays for each input.
For time-series models, nn = na, where na is the order of the AR model.
Note: The delay nk must be larger than 0. If you want nk = 0, shift the input sequence appropriately and use nk = 1 (see nkshift).
adm = 'ff' and adg = lam specify the forgetting factor algorithm with the forgetting factor λ=lam. This algorithm is also known as recursive least squares (RLS). In this case, the matrix P has the following interpretation: R2/2 * P is approximately equal to the covariance matrix of the estimated parameters.R2 is the variance of the innovations (the true prediction errors e(t)).
adm ='ug' and adg = gam specify the unnormalized gradient algorithm with gain gamma = gam. This algorithm is also known as the normalized least mean squares (LMS).
adm ='ng' and adg = gam specify the normalized gradient or normalized least mean squares (NLMS) algorithm. In these cases, P is not applicable.
adm ='kf' and adg =R1 specify the Kalman filter based algorithm with R2=1 and R1 = R1. If the variance of the innovations e(t) is not unity but R2; then R2* P is the covariance matrix of the parameter estimates, while R1 = R1 /R2 is the covariance matrix of the parameter changes.
Initial value of the parameters in a row vector, consistent with the rows of thm.
Default: All zeros.
Initial values of the scaled covariance matrix of the parameters.
Default: 104 times the identity matrix.
The argument phi0 contains the initial values of the data vector:
φ(t) = [y(t–1),...,y(t–na),u(t–1),...,u(t–nb–nk+1)]
If z = [y(1),u(1); ... ;y(N),u(N)], phi0 = φ(1) and phi = φ(N). For online use of rarx, use phi0, th0, and P0 as the previous outputs phi, thm (last row), and P.
Default: All zeros.
Estimated parameters of the model. The kth row of thm contains the parameters associated with time k; that is, the estimate parameters are based on the data in rows up to and including row k in z. Each row of thm contains the estimated parameters in the following order:
thm(k,:) = [a1,a2,...,ana,b1,...,bnb]
For a multiple-input model, the b are grouped by input. For example, the b parameters associated with the first input are listed first, and the b parameters associated with the second input are listed next.
Predicted value of the output, according to the current model; that is, row k of yhat contains the predicted value of y(k) based on all past data.
Adaptive noise canceling: The signal y contains a component that originates from a known signal r. Remove this component by recursively estimating the system that relates r to y using a sixth-order FIR model and the NLMS algorithm.
z = [y r]; [thm,noise] = rarx(z,[0 6 1],'ng',0.1); % noise is the adaptive estimate of the noise % component of y plot(y-noise)
If this is an online application, you can plot the best estimate of the signal y - noise at the same time as the data y and u become available, use the following code:
phi = zeros(6,1); P=1000*eye(6); th = zeros(1,6); axis([0 100 -2 2]); plot(0,0,'*'), hold on % Use a while loop while ~abort [y,r,abort] = readAD(time); [th,ns,P,phi] = rarx([y r],'ff',0.98,th,P,phi); plot(time,y-ns,'*') time = time + Dt end
This example uses a forgetting factor algorithm with a forgetting factor of 0.98. readAD is a function that reads the value of an A/D converter at the indicated time instant.