Identify continuous-time filter parameters from frequency response data

When building higher order models using high frequencies, it is important to scale the
frequencies, dividing by a factor such as half the highest frequency present in
`w`

, so as to obtain well-conditioned values of `a`

and
`b`

. This corresponds to a rescaling of time.

By default, `invfreqs`

uses an equation error method to identify the best
model from the data. This finds `b`

and `a`

in

$$\underset{b,a}{\mathrm{min}}{\displaystyle \sum _{k=1}^{n}wt(k){\left|h(k)A(w(k))-B(w(k))\right|}^{2}}$$

by creating a system of linear equations and solving them with the MATLAB^{®}
`\`

operator. Here
*A*(*w*(*k*)) and
*B*(*w*(*k*)) are the Fourier transforms
of the polynomials `a`

and `b`

, respectively, at the
frequency *w*(*k*), and *n* is the number
of frequency points (the length of `h`

and `w`

). This
algorithm is based on Levi [1]. Several variants have
been suggested in the literature, where the weighting function `wt`

gives
less attention to high frequencies.

The superior (“output-error”) algorithm uses the damped Gauss-Newton method for iterative search [2], with the output of the first algorithm as the initial estimate. This solves the direct problem of minimizing the weighted sum of the squared error between the actual and the desired frequency response points.

$$\underset{b,a}{\mathrm{min}}{\displaystyle \sum _{k=1}^{n}wt(k){\left|h(k)-\frac{B(w(k))}{A(w(k))}\right|}^{2}}$$

[1] Levi, E. C. “Complex-Curve
Fitting.” *IRE Trans. on Automatic Control*. Vol. AC-4, 1959,
pp. 37–44.

[2] Dennis, J. E., Jr., and R. B.
Schnabel. *Numerical Methods for Unconstrained Optimization and Nonlinear
Equations.* Englewood Cliffs, NJ: Prentice-Hall, 1983.