I misunderstood a little, "Invfreqs" tries to find a pole zero pattern that is close to the specified response but is stable itself.
But the conclusion still seems to be correct: that there is a direct calculation of the pole zero pattern, and the iterative additional calculation is just to approximate a more practical solution if the mathematically correct one is unrealisable.
Actually it calculates the numerator and denominator polynomial coefficients rather than the poles and zeros, but the raw coefficients are not very intuitive to me, so I write in terms of poles and zeros.
The example in the doco is an eye opener, there's no very obvious symmetry or transformation to connect the unstable solution and the approximation.
Which makes sense, if there was a simple transformation then we wouldn't need a iterative approximation solver.
David
Last edited:
To be clear, the original code was written in Fortran, back around 1985. That was translated to Matlab in 2006, and to C in 2008. In all cases, it was written in a very "explicit" fashion, to serve as a teaching tool. So there are no mathematical shortcuts, no fancy "modern programming" constructs, no optimizations.
Well, then, in spite of its age, it still feels more modern (vs. invfreqs), mainly because of its direct z-plane solution that appears to target the DSP-based implementations of the more recent times.
I wish that I could take credit for prescience, but honestly, at the time I just needed a thesis topic.
It's available directly from me. Just send a message to gjberchin (at) charter (dot) net, ask for FDLS Matlab code. There's also a C/C++ version, but I'm revising that now.
The past month issue of the AES journal has an open access review article by Balazs Bank on Warped, Kautz, and Fixed-Pole Parallel Filters (properties of invfreqz are also briefly discussed).
P.S.
Found closely related topic on Diyaudio
https://www.diyaudio.com/community/threads/python-open-room-correction-porc.215529/
https://www.aes.org/e-lib/browse.cfm?elib=21793
Frequency warping technique used to equalize the frequency response on an logarithmic scale looks impresive.P.S.
Found closely related topic on Diyaudio
https://www.diyaudio.com/community/threads/python-open-room-correction-porc.215529/
Last edited:
I agree: the mathematical technique is impressive. Using it produced an equalization curve that followed the peaks and dips of room reflections in the low end. Its's great that it could follow them but is filling in those nulls something we would want to do?
Perhaps I missed a mention earlier in the thread, but if a best-ft approximation is the aim it may be easier to go from the frequency response to its corresponding impulse response and use a Padé approximant to produce a rational function approximation. Dževad Belkić has published a number of papers and a book on related methods for analysing NMR data, REW uses that for extracting modal resonances. There is a summary of the method here.
In my experience, "go(ing) from the frequency response to its corresponding impulse response" is the hard part of that solution. If all that you have is a bunch of frequency response measurements, then how do you do that? In the continuous domain, if you have (or can find) the Laplace-domain transfer function corresponding to the measured frequency response, then the problem is already solved and there is no need for the Padé approximant. In the discrete domain, in order to compute the inverse Discrete Fourier (or similar) Transform you need frequency response measurements from DC to half the sampling frequency, equally-spaced at frequency steps sufficiently small (corresponding to a sufficiently long time record) to guarantee that there is no temporal aliasing. In practice, there is no guarantee that all of those criteria will be met.
From that page: "Frequency response is extrapolated to cover band from 0 Hz to Sample rate / 2 (Nyquist frequency) before inverse Fast Fourier Transform."
Exactly as I stated above. So they extrapolate to DC and half the sampling frequency, undoubtedly use interpolation to create the equal frequency response steps, and compute the inverse transform. That's a lot of approximation and assumption.
I did not read the code but the method, as posted by gberchin himself. But, I also see that that's already been answered.
Compared to what? The result will have the same frequency response as the original within the range of the original, isn't that the point? It's just a different representation of the same data.
It is up to the designer to decide how much and what kind of approximation error they are willing to tolerate.