Identify continuous-time filter parameters from frequency response data

`[b,a] = invfreqs(h,w,n,m)`

[b,a] = invfreqs(h,w,n,m,wt)

[b,a] = invfreqs(h,w,n,m,wt,iter)

[b,a] = invfreqs(h,w,n,m,wt,iter,tol)

[b,a] = invfreqs(h,w,n,m,wt,iter,tol,'trace')

[b,a] = invfreqs(h,w,'complex',n,m,...)

`invfreqs`

is the inverse operation of `freqs`

.
It finds a continuous-time transfer function that corresponds to a
given complex frequency response. From a laboratory analysis standpoint, `invfreqs`

is
useful in converting magnitude and phase data into transfer functions.

`[b,a] = invfreqs(h,w,n,m)`

returns
the real numerator and denominator coefficient vectors `b`

and `a`

of
the transfer function

$$H(s)=\frac{B(s)}{A(s)}=\frac{b(1){s}^{n}+b(2){s}^{n-1}+\cdots +b(n+1)}{a(1){s}^{m}+a(2){s}^{m-1}+\cdots +a(m+1)}$$

whose complex frequency response is given in vector `h`

at
the frequency points specified in vector `w`

. Scalars `n`

and `m `

specify
the desired orders of the numerator and denominator polynomials.

The length of `h`

must be the same as the length
of `w`

. `invfreqs`

uses `conj(h)`

at `-w`

to
ensure the proper frequency domain symmetry for a real filter.

`[b,a] = invfreqs(h,w,n,m,wt)`

weights
the fit-errors versus frequency, where `wt`

is a
vector of weighting factors the same length as `w`

.

`[b,a] = invfreqs(h,w,n,m,wt,iter)`

and

`[b,a] = invfreqs(h,w,n,m,wt,iter,tol)`

provide
a superior algorithm that guarantees stability of the resulting linear
system and searches for the best fit using a numerical, iterative
scheme. The `iter`

parameter tells `invfreqs`

to
end the iteration when the solution has converged, or after `iter`

iterations,
whichever comes first. `invfreqs`

defines convergence
as occurring when the norm of the (modified) gradient vector is less
than `tol`

, where `tol`

is an optional
parameter that defaults to 0.01. To obtain a weight vector of all
ones, use

invfreqs(h,w,n,m,[],iter,tol)

`[b,a] = invfreqs(h,w,n,m,wt,iter,tol,'trace')`

displays
a textual progress report of the iteration.

`[b,a] = invfreqs(h,w,'complex',n,m,...)`

creates
a complex filter. In this case no symmetry is enforced, and the frequency
is specified in radians between *–π* and *π*.

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.