System Identification Toolbox

Estimation from Frequency Domain Data

Here we give a very brief “cheat sheet” type explanation of how you estimate linear models to fit observed frequency domain data.

The toolbox supports frequency domain estimation with respect to both transfer function and state-space structures. In the transfer function case the model is either of output-error form:

$Y(omega_k) = G(gamma_k,theta) = frac{B(gamma_k)}{A(gamma_k)} + varepsilon_k$

or ARX form:

$A(gamma_k)Y(omega_k) = B(gamma_k) + varepsilon_k,qquad G(gamma_k,theta) = frac{B(gamma_k)}{A(gamma_k)}$

with varepsilon_k modelling measurement corruptions, and

$begin{array}{lcl} A(rho,theta) &=& 1+a_1rho^{-1}+a_2rho^{-2}+cdots+a_{n_a}rho^{-n_a} B(rho,theta) &=& b_0+b_1rho^{-1}+b_2rho^{-2}+cdots+b_{n_b}rho^{-n_b}. end{array}$

Here, rho = q and $gamma_k = e^{jomega_k}$ in the discrete time case, while rho = s and gamma_k = jomega_k in the continous time case.

In the state space case, the model is

$Y(omega_k) = G(gamma_k,theta) = C(gamma_{k}I-A)^{-1}B + D + varepsilon_k.$

Depending on which of these model structures is employed, The toolbox supports subpsace-based, least squares, and maximimum likelihood (ML) estimation methods.

First Fundamentals

Estimation of models via the command line involves executing the command

>> z = est(z,m,opt);

Here z specifies the data, m specifies the model structure, and opt determines optional aspects of how the estimate g is obtained. Because opt is optional, it is not strictly necessary to use it, so g=est(z,m) is perfectly legal.

Example data

In order to try out the above command using state-space model structures, we suggest you generate some test data. For example

fs=1;   N=1000;                                                              % Sampling frequency and number of measurements
den = real(poly([-0.1,-0.2,-0.02+j*1,-0.02-j*1,-0.01-j*0.1,-0.01+j*0.1]));   % Specify true system
num = 10*den(length(den));
w   = logspace(-3,log10(pi*fs),N);                                           % Choose measurement frequencies
F   = polyval(num,j*w)./polyval(den,j*w);                                    % Simulate the frequency response
Y   = F+1e-2*(randn(size(F)) + j*randn(size(F)));                            % Add some measurement noise to make interesting

Specifying the model structure

In the transfer function case, this is done by specifying the model orders n_a= m.A and (optionally) n_b= m.B (default is m.B=m.A). The toolbox then automatically assumes a transfer function structure is required, with defaulty m.type='oe’.

In the state space case, this is done by specifying the state dimension n_x= m.nx, which will then automatically specify a state space structure and set m.type='ss’.

Estimation Examples - Not Exhaustive

Discrete time ARX model Estimation via least squares
z.y = Y; z.w = w;          % Specify measurements and frequencies they were obtained at
m.A = 6;                   % Specify model order
m.type = 'arx';            % This is the default as well

g = est(z,m);              % Do the estimation
details(g);                % Present the results

Discrete time output-error model Estimation via least squares
z.y = Y; z.w = w;          % Specify measurements and frequencies they were obtained at
m.A = 6;                   % Specify model order
m.type = 'oe';             % Specify output error model

g = est(z,m);              % Do the estimation
details(g);                % Present the results

Continuous time output-error model Estimation via least squares
z.y = Y; z.w = w;          % Specify measurements and frequencies they were obtained at
m.A = 6;                   % Specify model order
m.type = 'oe';             % Specify output error model
m.op = 's';                % Specify a continuous time model

g = est(z,m);              % Do the estimation
details(g);                % Present the results

Continuous time state space model via Subspace Method
z.y = Y; z.w = w;          % Specify measurements and frequencies they were obtained at
m.nx = 6;                  % Specify state dimension
m.type = 'ss';             % Not strictly necessary - above implies it
m.op = 's';                % Specify a continuous time model
opt.alg = 'sid';           % Specify subspace method

g = est(z,m,opt);          % Do the estimation
details(g);                % Present the results

Continuous time state space model via Gauss-Newton search
z.y = Y; z.w = w;          % Specify measurements and frequencies they were obtained at
m.nx = 6;                  % Specify state dimension
m.type = 'ss';             % Not strictly necessary - above implies it
m.op = 's';                % Specify a continuous time model
opt.alg = 'gn';            % Specify Gauss--Newton Search

g = est(z,m,opt);          % Do the estimation
details(g);                % Present the results