An FFT-Based Equation-Error Method

The algorithm below minimizes the equation error in the frequency-domain.
As a result, it can make use of the FFT for speed. This algorithm is
implemented in Matlab's `invfreqz()` function when no iteration-count
is specified. (The iteration count gives that many iterations of the
Steiglitz-McBride algorithm, thus transforming equation error to output
error after a few iterations. There is also a time-domain implementation of
the Steiglitz-McBride algorithm called `stmcb()` in the Matlab Signal
Processing Toolbox, which takes the desired impulse response as input.)

Given a desired spectrum at equally spaced frequencies , with a power of , it is desired to find a rational digital filter with zeros and poles,

normalized by , such that

is minimized.

Since
is a quadratic form, the solution is readily obtained by
equating the gradient to zero. An easier derivation follows from minimizing
equation error variance in the time domain and making use of the orthogonality
principle [36].
This may be viewed as a system identification problem where the
known input signal is an impulse, and the known output is the desired
impulse response. A formulation employing an *arbitrary* known input
is valuable for introducing complex weighting across the frequency grid,
and this general form is presented. A detailed derivation appears in
[78, Chapter 2], and here only the final algorithm is given:

Given spectral output samples and input samples , we minimize

If is to be used as a weighting function in the filter-design problem, then we set .

Let : denote the column vector determined by , for filled in from top to bottom, and let : denote the size symmetric Toeplitz matrix consisting of : in its first column. A nonsymmetric Toeplitz matrix may be specified by its first column and row, and we use the notation :: to denote the by Toeplitz matrix with left-most column : and top row : . The inverse Fourier transform of is defined as

FFT

The scaling by is optional since it has no effect on the solution. We require three correlation functions involving and ,

where the overbar denotes complex conjugation, and four corresponding Toeplitz matrices,

where negative indices are to be interpreted mod
, *e.g.*,
.

The solution is then

where

[How to cite this work] [Order a printed hardcopy] [Comment on this page via email]

Copyright ©

Center for Computer Research in Music and Acoustics (CCRMA), Stanford University