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
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