matlab - Fitting an IIR filter to a complex transfer function

anybody knows of better ways to fit an IIR filter to a complex transfer function than Matlab's invfreqz.m? I'm implementing a little transfer function fitting toolbox for myself. I started with re-implementing the methods used in Matlab's invfreqz.m. Invfreqz.m is using Levy's complex curve fitting algorithm to compute starting coefficients for a iterative Gauss-Newton algorithm. Since Levy's paper is not documenting the algorithm very well and it's quite old, I'm looking for easier implementations and maybe even better algorithms to do the same job. Literature recherche turns out to be a real pain - so any hints from experienced people regarding literature, starting points or examples for implementations would be very appreciated.

Answer

As I mentioned in a comment, I think that the equation error method is a very good starting point for designing an IIR filter in the frequency domain with prescribed magnitude and phase responses. For an explanation of this method have a look at this answer and the references there.

Two problems associated with the equation error method are

1. stability of the designed filter is not guaranteed


2. the frequency domain error is weighted with the squared magnitude of the denominator polynomial of the designed filter, which means that the approximation is better in certain frequency regions than in others, but this weighting is a result of the design process and is usually undesirable.

The first problem can be tackled by adapting the total delay of the desired frequency response such that it matches with what a causal and stable filter with the chosen orders of numerator and denominator polynomials can achieve. So the desired complex frequency response $D(\omega)$ is changed to

$$\tilde{D}(\omega)=D(\omega)e^{-j\omega\tau}\tag{1}$$

with some delay parameter $\tau$ (which can also be negative). Optimizing the value of $\tau$ will help designing a stable filter with the smallest possible error. The only difference with the original desired response is an added (or subtracted) delay (i.e., a linear phase component).

One option to address the second problem is to solve a sequence of equation error problems with an adapted weight function, which is the basic idea of the Steiglitz-McBride method.

I've implemented the equation error method for the frequency domain design of IIR filters, and I've also added the possibility to iteratively approximate the optimal least squares solution, as described above. Note that convergence cannot be guaranteed, but in practice convergence is usually achieved after a few iterations. Also note that the resulting filter may be unstable, but this problem can be solved by adjusting the total delay of the desired response, as explained above. The Matlab/Octave code can be found at GitHub.

Here is a simple filter design example. The desired response is a low-pass filter with a linear phase response in the pass band. The delay is desired to be an integer number of samples, but the exact value is to be optimized to get the best approximation by a stable filter. The numerator and denominator orders are chosen to be $M=N=10$. Experiments showed that the best value for the desired pass band group delay is $28$ samples:

w = pi*[linspace(0,.2,20),linspace(.25,1,75)];   % frequency grid
D = [exp(-j*w(1:20)*28),zeros(1,75)];           % desired response
W = [ones(1,20),10*ones(1,75)];                  % error weighting
[a,b] = eqnerror(10,10,w,D,W,6);                 % 6 iterations of equation error method

The figure below shows the design result. Note that the pass band group delay approximates the desired $28$ samples. The designed filter is stable and has a maximum pole radius of $0.98$.