3

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.

lennon310
  • 3,590
  • 19
  • 24
  • 27
user967493
  • 387
  • 2
  • 12
  • 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. – Matt L. Feb 26 '16 at 14:00
  • 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? – user967493 Feb 26 '16 at 14:36
  • $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). – Matt L. Feb 26 '16 at 17:04
  • This answer is also relevant to your question. – Matt L. Feb 27 '16 at 09:47
  • 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... – user967493 Feb 27 '16 at 12:21
  • In my thesis you can find the (Matlab) code for all algorithms I discussed. What kind of filters are you interested in? – Matt L. Feb 27 '16 at 12:28
  • basically I want to fit a "relatively simple" transfer function with a IIR filter of a preferably low order. Thereby it's important to also approximate the phase response in the best way possible. As i used invfreqz.m in all my fittings before trying to write my own functions, i thought it was a good starting point, to get into the implementation. But as I already observed in your thesis, it doesnt turn out to be as good as i expected. – user967493 Feb 27 '16 at 12:36
  • This paper has a description of the Levy method (or one very similar to it) in the appendix: http://www.aes.org/e-lib/browse.cfm?elib=12429 – kippertoffee Feb 29 '16 at 09:02
  • thank you kippertoffee, that paper describes the alorithm very nicely – user967493 Feb 29 '16 at 09:19
  • Thank for your interesting discussion and the deep work you have done on the framework of your Ph.D; I have tried to use your IIR design filter function "mpiir_l2" on matlab which use "qp" function of the optimization toolbox. Unfortunately, this function is no longer part of this toolbox. I guess it has been relpaced by another function which is making the same job. Can you (or anybody else) let me know the name of it? Thanks very much in advance or how can i make your routine working? Thanks very much in advance. Warm regards Bernard – Bernard Jun 23 '21 at 15:02
  • @Bernard Welcome to SE.SP. Please do not post questions as answers. It is not the right thing to do on any Stack Exchange site. Google is your friend regarding your query.. – Peter K. Jun 23 '21 at 15:53

1 Answers1

3

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

Community
  • 1
Matt L.
  • 89,963
  • 9
  • 79
  • 179
  • 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? – user967493 Feb 28 '16 at 16:36
  • oh Matlab's mldivide() solves in Least squares sense. :) – user967493 Feb 28 '16 at 21:11
  • @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. – Matt L. Feb 29 '16 at 06:39
  • 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? – user967493 Feb 29 '16 at 09:43
  • @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. – Matt L. Feb 29 '16 at 10:19
  • aaaah now i see. – user967493 Feb 29 '16 at 10:26
  • Great answer and code, thanks! One thing I'm confused about (due to limited matlab and math knowledge): the least-squares solver in line 38, IIUC, you're solving for Ax = D, where A is the polynomial terms, with denominator pre-multiplied by the desired response (i.e. the unsquared RHS of E1 = |HkAk - Bk|^2, sans coefficients, in your other answer), D is your desired response, and x is the missing coefficients. However don't we want coefficients that bring E1 as close to zero as possible? I.e. shouldn't we be solving for Ax = E, where E is zeros, something like x = A\zeros()? – bryhoyt Apr 24 '17 at 23:35
  • And one other small thing I don't quite understand: why do you separate out the real and imaginary components when solving the equation? I.e. why x = [real(A);imag(A)]\[real(D);imag(D)] instead of x = A\D? – bryhoyt Apr 24 '17 at 23:36