elis, elisqa, elrpf2v, elrpv2f (Frequency Domain Identification Toolbox)

Basic iteration routine to calculate parametric estimate of linear transfer functions (elis); treat run parameters (elisqa, elrpf2v, elrpv2f).

Syntax

elis(Fdat,vdat,rppar)
[pvect,fit,Cp,CR,cfv] = ...
     elis(Fdat,vdat,rppar,fixp,rpalg,rppl,initp,rpfs)

elisqa
rpfout = elisqa(rpf,defaults)
elrpf2v(rpfile)
[rppar,fixp,rpalg,rppl,initp,rpfs] = elrpf2v(rpfile)
elrpv2f(rpfile,rppar)
elrpv2f(rpfile,rppar,fixp,rpalg,rppl,initp,rpfs)

Description

elis is the routine which performs the desired nonlinear least squares iteration to obtain the parameter estimates. elisqa is an optional question/answer routine. elis is rather complex, and the possibilities may generally need some explanation; in elisqa such explanations are given, and run parameters may be set. elrpf2v and elrpv2f perform conversions between run parameter vectors and a run parameter file.

The run parameters of elis can be set by its input arguments, or they can be passed through a so-called run parameter file, generated by elisqa. Most run parameters have default values, thus it is not necessary to give the values of all parameters for every run.

Fdat contains the input and output Fourier amplitudes: it may be defined as [freqv,x,y] (N-by-3 array), where freqv is the frequency vector, x is the complex input vector, and y is the complex output vector; or it may be a compound Fourier vector, generated by expfou; or it may be the name of a Fourier file.

vdat sets the variance values. It may be an N-by-2 or N-by-3 array, [varx,vary] or [varx,vary,covxy]; an 1-by-2 or 1-by-3 row vector, containing the constant variances; a variance vector (see expvar); or it may be a string with the name of a variance file.

Instead of either Fdat or vdat, the name of a run parameter file may also be given: this must be a string which ends by '.ebn'. It may not contain any other period, and the given run parameter file must define the Fourier or the variance data, respectively.

rppar, rpalg, rppl, rpfs are run parameter vectors, each composed of a set of run parameters. These prescribe the running of the algorithm, by elements that are either letters or numbers. It does not matter whether the vectors are given as strings or a numerical vectors, elis will take care of the proper meanings. The vectors may be shorter than their maximum length, or may contain NaN elements; in such a case the default values will be assigned to the corresponding run parameters. Some parameters may be influenced in different ways (e. g., the setting of starting values by rpalg(2) and by initp); in such cases the last setting will be valid, e. g., the one defined by initp in the above case.

rppar is the vector of the most often changed run parameters, associated with the model structure.

rppar(1) = domain: 's' or 'z'

rppar(2) = numord: order of the numerator

rppar(3) = denomord: order of the denominator

rppar(4) = fs: scaling angular frequency (in s-domain) or sampling frequency (in z-domain)

rppar(5) = 'a' for allpass design.

In the absence of a run parameter file, at least rppar(1:3) must be given.

fixp defines the fixed parameters. If it is given as an M-by-2 array, the first column must contain the serial numbers of the fixed parameters in the vector [num,denom,delay]', and the second column the values. For example, if the numerator order is 4 and the denominator order is 8, the parameter vector has (4+1+8+1+1=15) elements. The first denominator coefficient and the delay can be fixed by fixp = [6,1;15,0]. The delay alone can be fixed by fixp = [15,0].

If fixp is empty, previously set fixed parameters (by default or by a run parameter file) will not be changed.

A special form is when fixp is just a one-character string: 'n' for no fixed parameters at all (the delay is variable, too); 'f' for fixed delay, or 'v' for variable delay (for these two values the fixing of numerator and denominator coefficients remains untouched). Another special form is fixp = '0', for which the delay is fixed (to zero, or to the value defined by initp) and the zero-order coefficient of the denominator is set to 1.

rpalg provides control over the iteration possibilities.

rpalg(1) is the type of the iteration algorithm

'g' (Newton-Gauss),

'l' (Levenberg-Marquardt),

'm' (LM with svd),

's' (singular value decomposition),

'r' (Newton-Raphson) rpalg(2) = initset, way of setting the initial values:

'l' (ordinary least squares),

'w' (weighted least squares),

's' (WLS with svd),

'f' (file or parameter vector, see initp),

'e' (equation error method of the Signal Processing Toolbox) rpalg(3) = itmax: maximum number of iterations

rpalg(4) = rcostvar: stop if relative variation of cost function is smaller than this value.

rpalg(5) = rparvar: stop if relative variation of all parameters is smaller than this value.

Levenberg-Marquardt settings:

rpalg(6) = lambdadecr: number of consecutive decreases of the cost function before trial with lambda = 0

rpalg(7) = lambda: initial value of lambda

rpalg(8) = lambdalim: minimum value of lambda which allows stopping of iteration for small variations of the cost function or of the parameters.

rppl controls the plots during the iterations:

rppl(1) = plotdens: make plots after every plotdens cycle and after the initial setting and the last cycle if plotdens is finite; if it is inf, no plots will be made at all; if it is negative, the procedure is the same as for a positive number, except that the result of the initial setting will not be plotted.

rppl(2:5) = [fmin,fmax,amin,amax]: axis vector for the magnitude plots

rppl(6) = plotmode: type of frequency axis, 'i' for linear, 'o' for logarithmic

rppl(7) = calcrnum: calculate zeros from numerator polynomial, 'c' or 'n'

rppl(8) = calcrdenom: calculate poles from denominator polynomial, 'c' or 'n'

rppl(9:12): axis vector for zero/pole plots. Special meanings (the strings are always 4-character long, with appropriately set trailing spaces): [n,n,n,n] (four identical numbers): limit plotted zeros/poles by n*2*

*max(freqv) in the s-domain, or by just n in the z-domain; 'a' to show all zeros and poles; 'p' to show all zeros and poles, with the same scaling on the two axes.

rppl(13) = pzfollown: follow movement of zeros and poles by plotting several zero/pole sets, the number given by pzfollown.

initp may have different meanings: it may be the parameter vector of the initial parameter values; or a parameter filename; if it is a scalar number, it is the new value of the delay.

rpfs may contain the names of files to be generated, in the form of a string array, with trailing spaces if the filenames have different lengths. Possibilities: *.rep for elis a report file, *.par, *.pbn or *.pnt for generation of a parameter file, and *.cov, *.cbn or *.cnt for a covariance file. Other extensions are not allowed in rpfs.

If no argument is given for elis, elisqa will be started, with the name elisrpar.ebn.

elisqa offers default answers to the questions, these can be accepted by simply pressing Enter or Return.

The output arguments of elis are as follows.

pvect is the calculated parameter vector (see exppar).

fit is a vector containing information relevant to the results of the fit:

fit(1) = cost function
fit(2) = theoretical value of the cost function (c2/2)
fit(3) = 95% point of the theoretical distribution of the cf
fit(4) = number of frequencies
fit(5) = number of free parameters
fit(6) = performed iterations
fit(7) = last change of the cost function
fit(8) = last maximum of the relative changes of parameters
fit(9) = last lambda (for Levenberg-Marquardt, otherwise NaN)
fit(10) = approximate mean model error
fit(11) = mean absolute value of the transfer function
fit(12) = Akaike criterion
fit(13) = condition number of the matrix actually decomposed 
          or inverted in the last iteration step
fit(14) = condition number of J
fit(15) = condition number of Q=d^2C/dP^2, inverted when 
          calculating the approximate covariance matrix 
          (see Eq. (2.23)).
fit(16) = scaling frequency for internal calculations

The meaning and significance of the above condition numbers is explained in "Numerical Stability and Speed of the Procedures" in the User's Guide.

The approximate mean model error is calculated according to "Detection of Undermodeling and Overmodeling" in the User's Guide. The value of hmean is given, which will be imaginary if the cost function is too small (e. g., because the variances were overestimated, or simply because the cost function may be somewhat smaller than its expected value). This can be best checked from the information given in the report file. Hek and Xk are usually not even nearly constant. In these cases the square of the mean model error is averaged over all the frequencies given.

The mean model error can be compared to the mean absolute value of the transfer function, computed from the values at the same frequency points as above.

Cp is the approximate covariance matrix of the parameters (see "Covariance of the Estimate" in the User's Guide). CR is the approximate Cramér-Rao lower bound for the covariance matrix of the estimates. Cp and CR are calculated from the Jacobian of the last iteration (therefore at least one iteration step is necessary to calculate them). The covariances are usable only if the algorithm has converged.

cfv is the vector of the values of the cost function in each cycle (the initial cycle included). In the case of the Levenberg-Marquardt method cfv has a second column: the values of lambda are also given.

elis is a rather complex function, and the results often need careful documentation and studying. Therefore, a so-called report file can be requested with extensive textual information on the run (see rpfs above).

elisqa can be run separately, in order to set run parameters in a file for elis. The default run parameters are taken from the file defaults, or if this is not given, from the file rpf, or if neither of them is given, from an internal table (see below). If rpf is not given, the name of the file to be generated will be elisrpar.ebn.

The run parameter files are MATLAB binary files, with the extension '.ebn'. filenames given without extension will be extended in elisqa by '.ebn'. In principle, these run parameter files might be modified directly (e. g., by the routine savevar), but this is not recommended, since it is easy to make a mistake when doing this. elisqa offers a safe and easy way for such modifications.

The user of the routines usually need not bother about the internal names and values of the run parameters. However, the internal run parameters are listed below with their default values.

Sometimes it may be useful to use the same parameters in vector form, as given in a run parameter file, or vice versa. elrpf2v and elrpv2f serve this purpose. The meanings of the input and output arguments are explained above.

When elrpf2v has no output argument, or rpfile is empty when invoking elrpv2f, the values of the run parameters will be displayed on the screen.

The following paragraphs give a short description of the possibilities and solutions of elis, providing more detailed information than was possible in the description of the run parameter vectors.

The measurement data are supposed to be given by Fdat, or maybe in a Fourier file. For the internal calculations the frequency vector is scaled by the sampling frequency in the z-domain, or by a scaling frequency in the s-domain. The default setting for the scaling frequency is (

min +

max)/2.

When defining the model to be fitted, the domain must be specified (s or z), the orders of the numerator and the denominator must to be given, and the fixed parameters have to be defined (in the transfer function something must be fixed, since multiplication of each coefficient by the same constant gives the same transfer function). If no fixed nonzero parameters are given, elis will set the norm (the square root of the sum of squares) of the scaled coefficients (but not the delay) to 1. Another possibility is to set at least one coefficient of either the numerator or the denominator to a fixed value.

The allpass filter is treated as a special case of parameter fixing (the parameters of the denominator are the same as those of the numerator, but in reverse order), however, in this case at least one parameter of the numerator has to be fixed.

You can also specify whether the roots of the numerator and the denominator are to be calculated. When dealing with high order systems (>30), the iteration speed can be increased by sacrificing the pole/zero plot. However, the most significant acceleration can be achieved by completely suppressing plots (plotdens set to inf), but this has the risk of missing something that could be seen from the plot. A possible compromise is to set plotdens to a high value; in this case the starting values and the result of the last cycle will be plotted (independently of the actual serial number of the last cycle).

In the transfer function an extra delay term can be present. The value of this delay can be either fixed or estimated in the iteration procedure. However, a guess of the delay has to be given even if it is to be estimated, since the initial value may seriously influence the convergence properties.

The numerical methods solving the nonlinear LS problem (see "Basic Concepts" in the User's Guide) are standard methods in numerical analysis, including Newton-Gauss, Newton-Raphson, singular value decomposition (for the Newton-Gauss formulation), Levenberg-Marquardt and Levenberg-Marquardt with singular value decomposition. In the Levenberg-Marquardt algorithm an identity matrix, multiplied by lambda and by the Frobenius norm of JTJ, is added to JTJ before inversion; if a better fit is found, lambda is divided by 2, otherwise it is multiplied by 10 and the previous parameter values are restored. After lambdadecr successive divisions, an attempt is made with lambda = 0. (This corresponds to a Newton-Gauss step).

The singular value decomposition is always applied to J, even if it is combined with Levenberg-Marquardt, in order to make full use of the power of the singular value decomposition algorithm.

There are a few different possibilities to set the initial values for iteration. The first one is the standard initial setting procedure of ELiS: the linear least squares fitting. This can be modified to a weighted LS problem (weighting by the variances of the complex output amplitudes), and can also be solved by singular value decomposition. There are two more possibilities: the initial values can be taken from a parameter vector or file (e. g., the result of an earlier fit or of another method can be used), and also the equation error method of the Signal Processing Toolbox can be applied for initial value determination for elis (see invfreqz or invfreqs). In this latter case the input noise is transformed to the output as if a compound output noise was present.

For the control of the iterations, the maximum number of iteration cycles, the minimum relative variations of the cost function and of the estimated parameters can be set. The iteration will be terminated when the maximum number of iteration cycles is reached, or when any of the absolute values of the above variations is less than the corresponding minimum value. In the case of the Levenberg-Marquardt method the value of lambda is also considered: the iteration is terminated because of small variations only if the value of lambda is also smaller than lambdalim. The default value (1e-10) corresponds to 30 halvings of the default starting value 0.1.

A typical plot of the function elis is shown in the figure.

The left-hand side illustrates the fitting of the transfer function. The diagram consists of two parts. The upper part shows the absolute values of the ratios of the complex output and input amplitudes (+ marks), along with the magnitude of the estimated transfer function. The lower part shows the phase errors between the ratios of the measured complex amplitudes and the estimated transfer function. The plot is scaled vertically to (-180š,180š).

As a default, the pole/zero plot is also displayed at the right-hand side. The numbers of poles and zeros, the numbers of non-minimal phase zeros and unstable poles are also given, along with the number of zeros/poles occasionally not shown in the plot.

When elis is terminated, the current axes is the pole/zero plot. This can be manually rescaled if desired, but the text about unstable poles etc., has to be deleted first:

h = get(gca,'ch'); delete(h(1))

A somewhat simpler solution is to reapply plotelpz with a given axis vector:

plotelpz(pvect,axv,'','','nomsg')

If desired, even the uncertainty ellipses can be added:

plotelpz(pvect,axv,Cp,2,'nomsg')

The screen contains some further information concerning the run: the value of the cost function, the last change of the cost function, the theoretical expected value of the cost function (number of frequencies minus half of the number of free parameters), cycle number, iteration algorithm and way of setting the initial conditions, and the value of the delay.

A report file, that can be requested when running elis, will contain something like this:

FREQUENCY DOMAIN SYSTEM IDENTIFICATION TOOLBOX FOR MATLAB
ELiS RUN PROTOCOL, date and time: 1-Nov-93, 18:55:1
Report file: elisrpar.rep
Default run parameters, modified in command line

Fourier data given in command line
Experiment: 1, number of frequencies: 30
Input and output variances:
   5e-07   1.4754e-07
Input-output covariances are not given
Fit in s-domain, frequencies normalized internally
   by omegasc = 2.0452 rad*Hz
   suggestion: omegasc = 1.4408 rad*Hz
Orders: 1/3
No fixed nonzero parameters (norm=1 solution)
Fixed value of the delay:  0 s
Algorithm: Singular value decomposition
Initial value setting: WLS, by singular value decomposition

Allowed maximum number of iterations: 50,
   iterations performed: 4
Total run time: 0.56 min, time used for plots: 0.51 min
   time used for pole/zero calculations: 0.00 min
Stop if relative change of cost function is smaller than
   1.00e-06,  last relative variation: -1.59e-14
Stop if maximum relative change of parameters is smaller
   than 0.00e+00,  last max. rel. variation: +2.33e-11
Condition number of the decomposed or inverted matrix:
                       54.176, reciprocal: 0.018459
Condition number of J: 54.176, reciprocal: 0.018459
Condition number of Q=d^2C/dp^2:
                       2932.8, reciprocal: 0.00034097
                                       eps = 2.220e-16
Value of the cost function: 30.785, its double: 61.57
Theoretical value of the cost function: 27.0
Degrees of freedom of the chi-square value (2*cfth): 54
5%-95% points of the theoretical distribution of the cf:
   17.789, 38.098
Approximate mean model error: 0.00027292
Mean absolute value of the transfer function: 0.40298
Approximate relative mean model error: 0.00067726
Akaike criterion: 73.57
Number of free parameters: 6

Values of the parameters (with standard deviations, 
   calculated from the approximate covariance matrix):
Numerator:
s^0   -7.035487766480046e-02    std: 1.7580e-04  (0.24987%)
s^1   -7.038203864830019e-02    std: 5.6617e-05  (0.080443%)
Denominator:
s^0   -2.813969892740070e-01    std: 5.4478e-04  (0.1936%)
s^1   -2.110973116935604e-01    std: 1.6580e-04  (0.07854%)
s^2   -1.407034079272249e-01    std: 1.8950e-04  (0.13468%)
s^3   -7.037131510238287e-02    std: 9.1545e-05  (0.13009%)
Delay:
0 s    fixed

Zeros (rad*Hz):
-9.9961e-01

Poles (rad*Hz):
-1.6501e+00
-1.7467e-01 -1.5469e+00*j
-1.7467e-01 +1.5469e+00*j

First nonzero numerator coefficient is b(1)=-7.0382e-02
First nonzero denominator coefficient is a(3)=-7.0371e-02
b(1)/a(3)=1.0002e+00
Static gain: 2.5002e-01, -12 dB

Parameter file to save data: -
Covariance file to save data: -

%%%%% End of report file elisrpar.rep %%%%%

Default Argument Values

Here is the list of the run parameters and their default values (the possible answers are given between parentheses):

Ffile = ''; %name of Fourier file
vfile = ''; %name of variance file
pfile = ''; %parameter filename (maybe empty)
cfile = ''; %name of the covariance file to be generated
rfile = ''; %name of the report file to be generated
initpfile = ''; %name of param file, with initial values
algtype = 'NG'; %iteration algorithm (NG,LM,LMsvd,NR,svd)
amin = NaN; %axisvect(3) for magnitude plot
amax = NaN; %axisvect(4) for magnitude plot
calcrnum = 'c'; %calculate roots of numerator (c,n)
calcrdenom = 'c'; %calculate roots of denominator (c,n)
covxy = []; %input-output covariance
delay = 0; %(initial) value of the delay
delaytreat = 'f'; %delay fix (f) or variable (v)
denomord = 2; %order of the denominator
denomfix = []; %fixed denominator coefficients
denomfixind = []; %fixed denominator coefficient indices
domain = 's'; %domain of the model (s,z)
expi = 1; %serial number of the experiment to be used
fmin = 0; %lower bound of displayed frequencies
fmax = NaN; %max. displayed frequency in the s-domain
fs = NaN; %sampling (normalizing) frequency
initset = 'l'; %initial value setting (l,w,s,f,e)
itmax = 50; %maximum number of iteration cycles
lambda0 = .1; %starting lambda value for Levenberg-Marquardt
lambdadecr = 10; %after 10 consecutive decreases, 0 is tried
lambdalim = 1e-10; %iteration may stop below this value
numord = 1; %order of the numerator
numfix = []; %fixed numerator coefficients
numfixind = []; %fixed numerator coefficient indices
paramtreat = 'n'; %No fixed params (n), some are fixed:
   %(0, d, r) or allpass filter is designed (a)
plotdens = 1; %cycles of iteration to plot (1,2...inf)
pzfollown = 1; %pole/zero sets to be plotted on same plot
pzlimit = 'y'; %limit poles/zeros on plot to
   %pzlimitv*2**fmax (s) or 2 (z)
pzlimitv = 10; %limit poles/zeros on plot (see pzlimit)
rcostvar = 1e-6; %minimum relative variation of the cf
rparvar = 0; %minimum relative variation of parameters
vardef = 'u'; %variances: numbers or take from file (u,f)
varx = 1; %uniform variance value of input coefficients
vary = 1; %uniform variance value of output coefficients

Default values of the run parameter vectors:

rppar = [NaN,NaN,NaN,NaN];
rpalg = ['g','w',50,1e-6,0,10,0.1,1e-10];
rppl = [1,0,NaN,NaN,NaN,'i','c','c','a    ',1];
rpfs = '';

The NaN values mean that the actual values will be controlled by the Fourier data.

Examples

[freqv,x,y] = impfou('bandpass.fbn',1);
elis([freqv,x,y],3.4e-4*[1,1],['s',4,6]);
elis('inpchans.ebn');
elis('inpchanz.ebn',[],['z',16,16],[35,0]); %fixing the delay

elisqa('elisqtst.ebn','inpchans.ebn');

[rppar,fixp,rpalg,rppl,initp,prfs] = elrpf2v('inpchans.ebn');
elrpf2v('inpchans.ebn')

elrpv2f('newfile.ebn',['z',12,12],[0],'r');
elrpv2f('',['z',12,12],'0','r');

Diagnostics

The error and warning messages are self-explanatory. For the messages about condition numbers, see "Numerical Stability and Speed of the Procedures" in the User's Guide.

Algorithm

The algorithm and the main expressions are briefly described in Chapter 2, or in detail in [1].

References

[1] J. Schoukens and R. Pintelon, Identification of Linear Systems: a Practical Guideline for Accurate Modeling, London, Pergamon Press, 1991.

[ Previous | Help Desk | Next ]