0.13 Least squares design of iir filters  (Page 5/9)

Equation [link] represents an $N\times N$ system that can be solved for the coefficients $a_k$ given $c\left(n\right)$ for $0\le n\le N+M$ . These values can then be used in [link] to solve for the coefficients $b_k$ . The result is a system whose impulse response matches the first $N+M+1$ values of $f\left(n\right)$ .

Prony-based filter design methods

Both the original methods by Prony and Pade were meant to interpolate data from applications that have little in common with filter design. What is relevant to this work is their use of rational functions of polynomials as models for data, and the linearization process they both employ.

When designing FIR filters, a common approach is to take $L$ samples of the desired frequency response $D\left(\omega \right)$ and calculate the inverse DFT of the samples. This design approach is known as frequency sampling . It has been shown [link] that by designing a length- $L$ filter $h\left(n\right)$ via the frequency sampling method and symmetrically truncating $h\left(n\right)$ to $N$ values ( $N\ll L$ ) it is possible to obtain a least-squares optimal length- $N$ filter $h_N\left(n\right)$ . It is not possible however to extend completely this method to the IIR problem. This section presents an extension based on the methods by Prony and Pade, and illustrates the shortcomings of its application.

Consider the frequency response defined in [link] . One can choose $L$ equally spaced samples of $H\left(\omega \right)$ to obtain

where $A_k$ and $B_k$ represent the length- $L$ DFTs of the filter coefficients $a_n$ and $b_n$ respectively. The division in [link] is done point-by-point over the $L$ values of $A_k$ and $B_k$ . The objective is to use the relationship in described in [link] to calculate $a_n$ and $b_n$ .

One can express [link] as $B_k=H_k{A}_k$ . This operation represents the length- $L$ circular convolution $b\left(n\right)=h\left(n\right)◯\mathit{\text{L}}a\left(n\right)$ defined as follows [link]

where $h\left(n\right)$ is the length- $L$ inverse DFT of $H_k$ and the operator ${\left((·)\right)}_L$ represents modulo $L$ . Let

Therefore [link] can be posed as a matrix operation [link] of the form

where

is an $L\times L$ matrix. From [link] it is clear that the $L-(N+1)$ rightmost columns of $\widehat{\mathbf{H}}$ can be discarded (since the last $L-(N+1)$ values of $\widehat{a}$ in [link] are equal to 0). Therefore equation [link] can be rewritten as

or in matrix notation

where $a$ and $b$ correspond to the length- $N$ and $(M+1)$ filter coefficient vectors respectively and $\mathbf{H}$ contains the first $N+1$ columns of $\widehat{\mathbf{H}}$ . It is possible to uncouple the calculation of $a$ and $b$ from [link] by breaking $\mathbf{H}$ furthermore as follows,

Therefore

with

as defined in [link] . This formulation allows to uncouple the calculations for $a$ and $b$ using two systems,

Note that the last equation can be expressed as

where ${\mathbf{H}}_2=\left[{\widehat{h}}_2{\widehat{\mathbf{H}}}_2\right]$ (that is, ${\widehat{h}}_2$ and ${\widehat{\mathbf{H}}}_2$ contain the first and second through $N$ -th columns of ${\widehat{\mathbf{H}}}_2$ respectively).

From [link] one can conclude that if $L=N+M+1$ and if ${\widehat{\mathbf{H}}}_2$ and ${\mathbf{H}}_1$ are nonsingular, then they can be inverted In practice one should not invert the matrices ${\mathbf{H}}_1$ and ${\widehat{\mathbf{H}}}_2$ but use a more robust and efficient algorithm. See [link] for details. to solve for the filter coefficient vectors $a$ in [link] and solve for $b$ using ${\mathbf{H}}_1\tilde{a}=b$ .