Least-Squares Linear-Phase FIR Filter Design

Another versatile, effective, and often-used case is the weighted least squares method, which is implemented in the matlab function firls and others. A good general reference in this area is [204].

Let the FIR filter length be $L+1$ samples, with $L$ even, and suppose we'll initially design it to be centered about the time origin (``zero phase''). Then the frequency response is given on our frequency grid $\omega_k$ by

$$H(\omega_k) \eqsp \sum_{n=-L/2}^{L/2} h_n e^{-j\omega_kn}, \; k=0,1,2,\ldots,N-1, \; N\gg L.$$ (5.33)

Enforcing even symmetry in the impulse response, i.e., $h_n = h_{-n}$ , gives a zero-phase FIR filter that we can later right-shift $L/2$ samples to make a causal, linear phase filter. In this case, the frequency response reduces to a sum of cosines:

$$H( \omega_k ) \eqsp h_0 + 2\sum_{n=1}^{L/2} h_n \cos (\omega_k n), \quad k=0,1,2,\ldots, N-1,$$ (5.34)

or, in matrix form:

$$\underbrace{\left[ \begin{array}{c} H(\omega_0) \\ H(\omega_1) \\ \vdots \\ H(\omega_{N-1}) \end{array} \right]}_{{\underline{d}}} = \underbrace{\left[ \begin{array}{ccccc} 1 & 2\cos(\omega_0) & \dots & 2\cos[\omega_0(L/2)] \\ 1 & 2\cos(\omega_1) & \dots & 2\cos[\omega_1(L/2)] \\ \vdots & \vdots & & \vdots \\ 1 & 2\cos(\omega_{N-1}) & \dots & 2\cos[\omega_{N-1}(L/2)] \end{array} \right]}_\mathbf{A} \underbrace{\left[ \begin{array}{c} h_0 \\ h_1 \\ \vdots \\ h_{L/2} \end{array} \right]}_{{\underline{h}}} \protect$$ (5.35)

Recall from §3.13.8, that the Remez multiple exchange algorithm is based on this formulation internally. In that case, the left-hand-side includes the alternating error, and the frequency grid $\omega_k$ iteratively seeks the frequencies of maximum error--the so-called extremal frequencies.

In matrix notation, our filter-design problem can be stated as (cf. §3.13.8)

$$\min_{{\underline{h}}} \left\Vert \mathbf{A}{\underline{h}}-{\underline{d}}\right\Vert _2^2$$ (5.36)

where these quantities are defined in (4.35). We can denote the optimal least-squares solution by

$${\underline{\hat{h}}}\isdefs \arg \min_{\underline{h}}\left\Vert\,\mathbf{A}{\underline{h}}-{\underline{d}}\,\right\Vert _2 \eqsp \arg \min_{\underline{h}}\left\Vert\,\mathbf{A}{\underline{h}}-{\underline{d}}\,\right\Vert _2^2$$ (5.37)

To find ${\underline{\hat{h}}}$ , we need to minimize

$$\left\Vert\,\mathbf{A}{\underline{h}}-{\underline{d}}\,\right\Vert _2^2 = (\mathbf{A}{\underline{h}}-{\underline{d}})^T(\mathbf{A}{\underline{h}}-{\underline{d}})$$

$$= ({\underline{h}}^T\mathbf{A}^T-{\underline{d}}^T)(\mathbf{A}{\underline{h}}-{\underline{d}})$$

$$= {\underline{h}}^T\mathbf{A}^T\mathbf{A}{\underline{h}} -{\underline{h}}^T\mathbf{A}^T{\underline{d}} -{\underline{d}}^T\mathbf{A}{\underline{h}} +{\underline{d}}^T{\underline{d}} \protect$$ (5.38)

This is a quadratic form in ${\underline{h}}$ . Therefore, it has a global minimum which we can find by setting the gradient to zero, and solving for ${\underline{h}}$ .^5.14Assuming all quantities are real, equating the gradient to zero yields the so-called normal equations

$$\mathbf{A}^T\mathbf{A}{\underline{h}}\eqsp \mathbf{A}^T{\underline{d}}$$ (5.39)

with solution

$$\zbox {{\underline{\hat{h}}}\eqsp \left[(\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\right]{\underline{d}}.}$$ (5.40)

The matrix

$$\mathbf{A}^\dagger \isdefs (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T$$ (5.41)

is known as the (Moore-Penrose) pseudo-inverse of the matrix $\mathbf {A}$ . It can be interpreted as an orthogonal projection matrix, projecting ${\underline{d}}$ onto the column-space of $\mathbf {A}$ [264], as we illustrate further in the next section.