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.
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
- stability of the designed filter is not guaranteed
- 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$.
-
$\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 specifyom
. 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 useom
for $\omega$ ('omega'), and sometimes I usew
, and this time they got mixed up. I've corrected the code. $\endgroup$ – Matt L. Feb 29 '16 at 10:19
invfreqz.m
... $\endgroup$ – user967493 Feb 27 '16 at 12:21