Design $ {L}_{2} $ Norm Optimal Infinite Impulse Response (IIR) Filters

It is widely known that matching a FIR filter of fixed length to a band model is an unconstrained QP-problem. The MATLAB function firls() implements a solution to this problem. Basically, one symbolically evaluates integrals of the (Z-domain) transfer function over the unit circle, getting astonishingly simple expressions for the match quality. Additionally, that method is stable and will not fail, other than the Parks McClellan Filter Design Algorithm.

Is it possible to fit IIR filters to such a model in the $ {L}_{2} $ Norm Sense?
Several problems arise:

1. While the cost function (fit quality?) can be written down, is it not a simple quadratic anymore, but a fraction of two quadractics.


2. You appear to have to optimize nominator and denominator at the same time.


3. It might be hard to constrain the denominator to "minimum phase", so that the resulting filter is stable.

I am not looking for gradient descent or other local optimization methods, as the cost function might have several local optima. An idea that crossed my mind would be that it appears possible to pose the FIR problem so that it yields the optimal coefficients for an IIR all-pass.

Answer

As you have already pointed out in your question, it is not possible (without using optimization methods) to compute an exact L2 solution for the frequency domain design problem of IIR filters due to the non-linear relationship between the filter coefficients and the error function. There is, however, a method which can come close and which transforms the problem to a linear one: the Equation Error Method. Instead of defining the error measure as

$$E_0=\sum_k\left |H_k -\frac{B_k}{A_k}\right|^2\tag{1}$$

(where $k$ is the frequency index, $H_k$ is the complex desired frequency response at frequency $\omega_k$, and $B_k$ and $A_k$ are the numerator and denominator polynomials, respectively, also evaluated at the frequency point with index $k$), one defines an error measure

$$E_1=\sum_k\left |H_kA_k -B_k\right|^2\tag{2}$$

Minimizing $E_1$ is a linear problem in the filter coefficients. You get an overdetermined system of linear equations which can be solved (in the $l_2$ sense) by solving a set of linear equations. (Of course, the number of frequency points must be greater than the number of filter coefficients). The result of minimizing (2) is identical to solving a weighted $l_2$ problem with weight function $|A_k|^2$:

$$E_1=\sum_kW_k\left|H_k -\frac{B_k}{A_k}\right|^2\quad\textrm{with}\quad W_k=|A_k|^2$$

There are two (related) problems with this method:

1. 
   

   The specification must not only include the desired magnitude response but also the desired phase. If the phase specification is chosen in a way that does not fit the chosen filter order and the general properties of IIR filters, the approximation error will be large, and the filter might be unstable (which brings us to the next point).

   
   



2. 
   

   Stability is not considered in the design process. Depending on the specification, the best approximation could either be stable or unstable.

   
   

Of course, poles outside the unit circle can be reflected inside the unit circle without affecting the designed magnitude response, but the approximation error might still be large, because an unstable filter indicates that the specification is not well suited to the chosen filter order. So it takes quite some experience to choose an appropriate desired phase response.

An in-depth explanation of the equation error method can be found in the book Digital Filter Design by Parks and Burrus. You can find a good overview here.