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 $H(e^{j\omega_k})$ at equally spaced frequencies $\omega_k = 2\pi k/N, k=0,\ldots\,,N-1$ , with a power of , it is desired to find a rational digital filter with ${{n}_b}$ zeros and ${{n}_a}$ poles,

$\displaystyle \hat{H}(z) \isdef \frac{\hat{B}(z)}{\hat{A}(z)} \isdef \frac{\sum_{k=0}^{{n}_b}b_k z^{-k}}{\sum_{k=0}^{{n}_a}a_k z^{-k} },$

normalized by

, such that

$\displaystyle J^2_E = \sum_{k=0}^{N-1} \left\vert\hat{A}(e^{j\omega_k})H(e^{j\omega_k})-\hat{B}(e^{j\omega_k})\right\vert^2$

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 $Y(e^{j\omega_k})$ and input samples $U(e^{j\omega_k})$ , we minimize

$\displaystyle J^2_E = \sum_{k=0}^{N-1} \left\vert\hat{A}(e^{j\omega_k})Y(e^{j\omega_k})-\hat{B}(e^{j\omega_k})U(e^{j\omega_k})\right\vert^2 .$

If $\vert U(e^{j\omega_k})\vert^2$ is to be used as a weighting function in the filter-design problem, then we set $Y(e^{j\omega_k}) = H(e^{j\omega_k})U(e^{j\omega_k})$ .

Let $\underline{x}[n_1\,$ : $\,n_2]$ denote the column vector determined by , for $n=n_1,\ldots\,,n_2$ filled in from top to bottom, and let $T(x[n_1\,$ : $\,n_2])$ denote the size symmetric Toeplitz matrix consisting of $\underline{x}[n_1\,$ : $\,n_2]$ in its first column. A nonsymmetric Toeplitz matrix may be specified by its first column and row, and we use the notation $T(\underline{x}[n_1\,$ : $\,n_2],\underline{y}^T[m_1\,$ : $\,m_2])$ to denote the by Toeplitz matrix with left-most column $\underline{x}[n_1\,$ : $\,n_2]$ and top row $\underline{y}^T[m_1\,$ : $\,m_2]$ . The inverse Fourier transform of $X(e^{j\omega_k})$ is defined as

$\displaystyle x(n) =$ FFT $\displaystyle ^{-1}{X(e^{j\omega_k})} \isdef \frac{1}{N}\sum_{k=0}^{N-1} X(e^{j\omega_k}) e^{j\omega_k n} .$

The scaling by

is optional since it has no effect on the solution. We require three correlation functions involving

and

$\begin{eqnarray*} \underline{R}_{uu}(n) &\isdef & \mbox{FFT}^{-1}{\left\vert U(e... ...})\overline{U(e^{j\omega_k})}} \\ n & = & 0,1,\ldots\,,N-1,\\ \end{eqnarray*}$

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

$\begin{eqnarray*} R_{yy} &\isdef & T(\underline{R}_{yy}[0\,\mbox{:}\,{{n}_a}-1])... ...u}^T[-1\,\mbox{:}\,-{{n}_a}])\\ R_{uy} &\isdef & R_{uy}^T , \\ \end{eqnarray*}$

where negative indices are to be interpreted mod , e.g., $R_{yu}(-1)=R_{yu}(N-1)$ .

The solution is then

$\displaystyle \hat{\theta}^\ast = \left[\begin{array}{c} \underline{\hat{B}}^\a... ...,{{n}_b}] \\ [2pt] \underline{R}_{yy}[1\,\mbox{:}\,{{n}_a}] \end{array}\right]$

where

$\displaystyle \underline{\hat{B}}^\ast \isdef \left[\begin{array}{c} \hat{b}^\a... ...t{a}^\ast _0 \\ [2pt] \vdots \\ [2pt] \hat{a}^\ast _{{n}_a}\end{array}\right],$

Previous: Stability of Equation Error Designs
Next: Prony's Method

Order a Hardcopy of Introduction to Digital Filters

About the Author: Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.

Sign in

Search Online Books

Free Online Books

Ads

Chapters

See Also

Introduction to Digital Filters
    Recursive Digital Filter Design
       Filter Design by Minimizing the L2 Equation-Error Norm
          An FFT-Based Equation-Error Method