3.8 Direct frequency domain iir filter design methods  (Page 6/6)

Prony's method for time-domain design of iir filters

The problem of designing an IIR digital filter with a prescribed time-domain response is addressed in this section.Most formulations of time-domain design of IIR filters result in nonlinear equations for the same reasons as for frequency-domaindesign. Prony, in 1790, derived a special formulation for the analysis of elastic properties of gases, which resulted in linearequations. A more general form of Prony's method can be applied to the IIR filter design by use of a matrix description [link] , [link] .

The transfer function of an IIR filter is given by

and the impulse response $h\left(n\right)$ is related to $H\left(z\right)$ by the $z$ transform.

[link] can be written

which is the z-transform version of convolution. This convolution can be written as a matrix multiplication. Using the first K+1terms of the impulse response, this is written

In order to uncouple the calculations of the $a_n$ and the $b_n$ , the matrices are partitioned to give

where $\mathbf{b}$ is the vector of the $M+1$ numerator coefficients of [link] , ${\mathbf{a}}^{*}$ is the vector of the $N$ denominator coefficients ( $ao=1$ ), ${\mathbf{h}}_1$ is the vector of the last $(K-M)$ terms of the impulse response, ${\mathbf{H}}_{\mathbf{1}}$ is the $M+1$ by $N+1$ partition of [link] , and ${\mathbf{H}}_2$ is the $(K-M)$ by $N$ remaining part. The lower $K-M$ equations are written

or

which must be solved for ${\mathbf{a}}^{*}$ . The upper $M+1$ equations of [link] are written

which allows the calculation of $\mathbf{b}$ .

If $L=N+M$ , then ${\mathbf{H}}_2$ is square. If ${\mathbf{H}}_2$ is nonsingular, [link] can be solved exactly for the denominator coefficients in ${\mathbf{a}}^{*}$ , which are augmented by the unity term to give $\mathbf{a}$ . From [link] , the numerator coefficients in $\mathbf{b}$ are found. If ${\mathbf{H}}_2$ is singular [link] and there are multiple solutions, a lower order problem can be posed. If there are nosolutions, the methods of the next section must be used.

Note that any order numerator and denominator can be prescribed. If the filter is in fact an FIR filter, $\mathbf{a}$ is unity and ${\mathbf{a}}^{*}$ does not exist. Under these conditions, [link] states that $b_n=h_n$ , which is one of the cases of FIR frequency sampling covered in Section 3.1 of [link] . Also note that there is no control over the stability of the filter designed by thismethod.

Although Prony's method, applied to the time-domain design problem here, is similar to the solution of the frequency-samplingIIR design problem, there are important differences. The inverse DFT is used to obtain the matrix in the frequency domain problem, which iscyclic convolution. Equation [link] is noncyclic convolution and the $K+1$ terms of $h\left(n\right)$ , used to form $\mathbf{H}$ , result from a truncation of the infinitely long sequence.

An approximate solution or the least equation error problem

In order to obtain better practical filter designs, the interpolation scheme of the previous section is extended to give an approximation designmethod [link] . It should be noted at the outset that the method developed in this section minimizes an equation-error measure and not the usualfrequency-response error measure.

The number of samples specified, $L+1$ , will be made larger than the number of filter coefficients, $M+N+1$ . This means that ${\mathbf{H}}_2$ is rectangular and, therefore, [link] cannot in general be satisfied. To formulate an approximation problem, a length- $(L+1)$ error vector $\epsilon$ is introduced in [link] and [link] to give

[link] becomes

where now ${\mathbf{H}}_2$ is rectangular with $(L-M)>N$ . Using the same methods as used to derive [link] , the error $\epsilon$ is minimized in a least-squared error sense by the solution of the normal equations [link]

If the equations are not singular, the solution is

If the normal equations are singular, the pseudo-inverse [link] can be used to obtain a minimum norm or reduced order solution.

The numerator coefficients are found by the same techniques as before in [link]

which results in the upper $M+1$ terms in $\epsilon$ being zero and the total squared equation error being minimum.

As is true for LS-error design of FIR filters, [link] is often numerically ill-conditioned and [link] should not be used to solve for ${\mathbf{a}}^{*}$ . Special algorithms such as those used by Matlab and LINPACK [link] , [link] should be employed.

Various modifications can be made to the form of Prony's method presented. After the denominator is found by minimizingthe equation error, the numerator can be found by minimizing the solution error. It is possible to mix the exact and approximatemethods. The details can be found in [link] , [link] , [link] .

Several modifications to Prony's method have been made to use it to minimize the solution error. Most of these iterativelyminimize a weighted-equation error with Prony's method and update the weights from the previous determination of $a$ [link] , [link] .

If an LS-error, time-domain approximation is the desired result, a minimization technique can be applied directly to thesolution error. The most successful method seems to be the Gauss- Newton algorithm with a step-size control. This combined withProny's method to find starting parameters is an effective design tool.