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 $N$ a power of $2$ , it is desired to find a rational digital filter with ${{n}_b}$ zeros and ${{n}_a}$ poles,

$$\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 $a_0 = 1$ , such that

$$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 $J^2_E$ 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. More detailed background 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

$$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 $x(n)$ , for $n=n_1,\ldots\,,n_2$ filled in from top to bottom, and let $T(x[n_1\,$:$\,n_2])$ denote the size $n_2-n_1+1$ 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 $n_2-n_1+1$ by $m_2-m_1+1$ 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

$$x(n) =$$   FFT$$^{-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 $1/N$ is optional since it has no effect on the solution. We require three correlation functions involving $U$ and $Y$ ,

$$\begin{eqnarray*} \underline{r}_{uu}(n) &\isdef & \mbox{FFT}^{-1}{\left\vert U(e^{j\omega_k})\right\vert^2} \\ \underline{r}_{yy}(n) &\isdef & \mbox{FFT}^{-1}{\left\vert Y(e^{j\omega_k})\right\vert^2} \\ \underline{r}_{yu}(n) &\isdef & \mbox{FFT}^{-1}{Y(e^{j\omega_k})\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])\\ R_{uu} &\isdef & T(\underline{r}_{uu}[0\,\mbox{:}\,{{n}_b}])\\ R_{yu} &\isdef & T(\underline{r}_{yu}[-1\,\mbox{:}\,{{n}_b}-1],\underline{r}_{yu}^T[-1\,\mbox{:}\,-{{n}_a}])\\ R_{uy} &\isdef & R_{yu}^T , \\ \end{eqnarray*}$$

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

The solution is then

$$\hat{\theta}^\ast = \left[\begin{array}{c} \underline{\hat{B}}^\ast \\ [2pt] -\underline{\hat{A}}^\ast \end{array}\right] = \left[\begin{array}{cc} R_{uu} & R_{yu} \\ [2pt] R_{uy} & R_{yy} \end{array}\right]^{-1} \left[\begin{array}{c} \underline{r}_{yu}[0\,\mbox{:}\,{{n}_b}] \\ [2pt] \underline{r}_{yy}[1\,\mbox{:}\,{{n}_a}] \end{array}\right],$$

where

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