2
$\begingroup$

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.

Improve this question
$\endgroup$
  • $\begingroup$ Have a look at chapter 5 of this thesis and the references therein. Also have a look at this answer explaining the equation error method, which is one of the most basic methods to approximate a complex frequency response by an IIR filter. $\endgroup$ – Matt L. Feb 26 '16 at 14:00
  • $\begingroup$ thank you for your thesis! :) Seems like some great input, I will have a look into it and try to comprehend what you have done . And thank you aswell for the great explanation of the error equation method. Since I have already looked into the method, I have one remark to your explanation: the error equation method is based on the assumption that d(n) ≈ y(n) . So the E1 should differ slightly from the optimum E0. Or am I misunderstanding something? $\endgroup$ – user967493 Feb 26 '16 at 14:36
  • $\begingroup$ $E_1$ is a weighted version of $E_0$, where the weighting is equal to the squared magnitude of the denominator (as explained in the other answer). The version I explained in that answer is a frequency domain method, so there is no reference to $d[n]$ or $y[n]$ (which I interpret as time domain sequences). $\endgroup$ – Matt L. Feb 26 '16 at 17:04
  • $\begingroup$ This answer is also relevant to your question. $\endgroup$ – Matt L. Feb 27 '16 at 9:47
  • $\begingroup$ thank you for the great input Matt! Do you by coincidence know if there is some literature or web tutorials or anything available describing a more pracitcal approach? Reason is, the theory is to some extent clear after reading, but the actual implementation seems to be miles away from the theoretical approach. If you compare for example Levy's description of the algorithm to the acutal implementation in invfreqz.m... $\endgroup$ – user967493 Feb 27 '16 at 12:21
1
$\begingroup$

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$.

enter image description here

Improve this answer
$\endgroup$
  • $\begingroup$ first of all thank you a lot. i have some questions regarding the code: how do you solve the IIR design in an LMS sense? The following line in your code: 'x = [real(A);imag(A)][real(D);imag(D)];' is a very logic step. Most literature i have read, solve the LMS problem through usage of the autocorrelation and cross corelation. I suppose you minimize the error function through changing the weighting function with respect to the a coefficients in this line: 'W = W0./abs(den);' right? But how does it solve the error minimzation in the LMS sense? $\endgroup$ – user967493 Feb 28 '16 at 16:36
  • $\begingroup$ oh Matlab's mldivide() solves in Least squares sense. :) $\endgroup$ – user967493 Feb 28 '16 at 21:11
  • $\begingroup$ @user967493: The backslash operator '\' does the job of solving an overdetermined system of linear equations in the least squares sense. It is indeed equivalent to mldivide. $\endgroup$ – Matt L. Feb 29 '16 at 6:39
  • $\begingroup$ Again, thank you very much Matt L.! Just out of interest in your example, you set the desired function D to D = [exp(-j*om(1:20)*28),zeros(1,75)];. But you didn't specify om. What did you set it to? $\endgroup$ – user967493 Feb 29 '16 at 9:43
  • $\begingroup$ @user967493: Good catch, that's a mistake. It should be w. Sometimes I use om for $\omega$ ('omega'), and sometimes I use w, and this time they got mixed up. I've corrected the code. $\endgroup$ – Matt L. Feb 29 '16 at 10:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.