No, no, I meant the actual source-code for this function. Never mind, I am able to see the code from the following path:
...\MATLAB\<version>\toolbox\signal\signal\invfreqs.m
I'm attaching the file here with a TXT entension, for anyone interested in seeing it. Open with Wordpad to see the source code.
...\MATLAB\<version>\toolbox\signal\signal\invfreqs.m
I'm attaching the file here with a TXT entension, for anyone interested in seeing it. Open with Wordpad to see the source code.
My experience with invfreqs() has been hit-or-miss. Still, it's a viable candidate, worthy of at least a try.
gberchin:
Sir, in my opinion, there are two things to consider:
1) The MATLAB source code needs to be "knocked down" before actually understanding / implementing it. However, with your method, the OP has the inventor himself to guide him through the process of implementing the method. It is also time-tested and mathematically documented and since it has also been published, I assume it was peer reviewed as well.
2) I was able to see Mathworks' copyright somewhere within invfreqs().m, so I guess it will have to be re-written by the OP, before use. However, you've offered to unconditionally share your method (and related documentation) for the OP's benefit, which I believe makes things much much easier (copyright-wise) when compared to the Mathworks route.
Nevertheless, the usage of the continuous frequency s-domain is one solid advantage that the Mathworks' method might have. I think we'll wait and see what the OP decides.
Sir, in my opinion, there are two things to consider:
1) The MATLAB source code needs to be "knocked down" before actually understanding / implementing it. However, with your method, the OP has the inventor himself to guide him through the process of implementing the method. It is also time-tested and mathematically documented and since it has also been published, I assume it was peer reviewed as well.
2) I was able to see Mathworks' copyright somewhere within invfreqs().m, so I guess it will have to be re-written by the OP, before use. However, you've offered to unconditionally share your method (and related documentation) for the OP's benefit, which I believe makes things much much easier (copyright-wise) when compared to the Mathworks route.
Nevertheless, the usage of the continuous frequency s-domain is one solid advantage that the Mathworks' method might have. I think we'll wait and see what the OP decides.
1) The Matlab help page references two documents that they used to develop the invfreqs() function. I would recommend that one refer to those, rather than reverse-engineering the Matlab function.
2) FDLS has been freely available for decades, but it's not well-known. I help however I can, but it's not something to which I devote a lot of time any more.
As for s-domain vs. z-domain; if the equalization is to be done in DSP, then starting with a z-domain transfer function will probably make the job easier.
Oh, yes ! One could just directly move to formatting and copying the coefficients, bypassing bi-linear transform, warping compensation etc.
I purposely did not use the term "reverse-engineering", as I never meant hackers' disassembly of the function, but only an analysis to help one understand (learn) what is done. Of course, the mathematics behind the function is more easily understood from the citations.
Besides, I also believe the programmers at Mathworks to be more skilled than many of us and therefore consider it to be a good idea to visit the source code (to see if any extra "tricks" have been used), in addition to reading the mathematical explanation in the citations.
Comprehensive review by
Vesa Välimäki and Joshua D. Reiss, All About Audio Equalization: Solutions and Frontiers, Appl. Sci. 2016, 6(5), 129.
Also look for German Ramos papers in AES journal/conference proceedings on the subject.
Do you have an explanation why the Laplace version is worse?
That doesn't surprise me, kind of to be expected .
Thanks, I was unaware of this...
and unaware there was an explanation and references, I am more interested in the explanation than a particular solution, at the moment.
What problems have you seen? Do they illuminate any issues?
Thanks, mostly review but there are some useful details.
One of their references is the same as in the Mathworks "Invfreq", a paper by Levi. that looks to be the classic for the subject.
But it's in a journal from 1959 (!) and not easy to find. After 60 years I had hoped for public domain.
Best wishes
David
No. Once I found that it did not solve my problem, I moved on.
Exactly what one would expect -- sometimes it worked well, sometimes it didn't. It was not my intent to analyze the algorithm, but to solve my problem. So again, once I found that it did not, I moved on.
Another candidate: Oppenheim & Schafer mention a technique by Steiglitz in section 5.3.1 of their "Digital Signal Processing" book.
K. Steiglitz, "Computer-Aided Design of Recursive Digital Filters", IEEE Transactions on Audio and Electroacoustics, Vol. AU-18, June 1970.
I've never tried this method.
K. Steiglitz, "Computer-Aided Design of Recursive Digital Filters", IEEE Transactions on Audio and Electroacoustics, Vol. AU-18, June 1970.
I've never tried this method.
Last edited:
I glanced over the entire discussion, as far as I understood, you are generally interested in the problem of a complex function approximation (nonlinear real approximation problem). Along with the mentioned by gberchin Steiglitz method, there also is Chebyshev method for complex function approximation problem. You can find its description in T. W. Parks and C. S. Burrus, Digital Filter Design, John Wiley & Sons, New York, 1987 (chapter 4) (Fortran code for complex approximation is presented).
Last edited:
Thank you.
A cursory reading seems to indicate that Levi only uses the magnitude error in the algorithm. If that's the case, then the model may not match the phase response of the prototype system.
Maybe, but the MATLAB function combines the methods in both (cited) papers and seems to use Levi's method only to arrive at an educated guess for the next one, which I think, is what it finally implements.
Dave Zan:
Please note that the MATLAB function is iterative, assuming you're looking for direct, non-iterative methods.
Dave Zan:
Please note that the MATLAB function is iterative, assuming you're looking for direct, non-iterative methods.
Not problematic at all, provided some knowledge is held by the user.
I.e. the simple rule that 0.707 Q is flat, 0dB gain, in all cases.
Remember the simple rule and its easy
Sorry, that was supposed to be "shelf slope of 1" that got messed up during editing. But how about "a simple rule" for a gain & slope of 2 (each) ? I think that would really help. 😉
The "problem" lies in that some processing units ask the user for the Q-factor while some others the slope, and the relation between these quantities is just too "complicated" for most people to remember:
And, for anyone using RBJ's definition of shelving, A = sqrt(A) and the equation becomes:
Neither seems to be friendly from my point of view. I recently wrote some programmes for my "DIY Dolby Box", where I decided to ditch the quantity called slope and use Q-factor instead.
I'd agree with that.
Slope factor is just a term that confuses the situation.
A gain of 2, 6dB, should be straightforward to figure. Shouldn't it?
As for automating the iteration, that's a different story perhaps, and not my forte.
Mind you, Boxsim does this auto correction; albeit with a different proviso - minimise deviation from a flat acoustic summing of drivers - so curve fitting, to an extent, it's just that the curve, can't be changed.
Yes, I am a data analyst and statistician and this overlaps with my audio hobby.
In statistics the data is always real, and I have never seen any work on how to fit complex data.
So it interests me as a maths problem as well as its audio applications.
Thanks, I will try to track this down.
I think he includes the phase error in the error term and then tries to minimise that. But I could be mistaken, it's a bit cryptic.
I think you have perhaps misinterpreted here.
It seems the Matlab function calculates the poles and zeroes of the function directly and that allows calculation of the inverse.
But if the function has non minimum phase (RHP zeroes) then the inverse will not be stable.
That's just an inconvenient reality.
So the iterative step is only required in this case, to try to find an approximation such that the inverse remains stable.
As you say, it calculates Levy's solution curve fit then uses it as the start of the iterative approximation, but this doesn't invalidate the fact that the curve fit can be directly calculated.
I could, of course, be mistaken here too.
Best wishes
David
Last edited:
Quite possibly, but the following lines seem to suggest that the iterative procedure is executed whenever more than three arguments (besides h and w) are entered by the user, for example:
invfreqs(h, w, n, m, [], iter, tol);
Code:
gaussFlag = length(varargin)>3; % run Gauss-Newton algorithm or not?
......
if ~gaussFlag,return,end
% Now for the iterative minimization
if isempty(maxiter), maxiter = 30; end
if isempty(tol)
tol=0.01;
end
.....If the idea behind all this is equalisation, then it may be better to deal with zeroes in the following manner and avoid instability, as we've seen earlier in this thread. The poles (being on the LHS) may either be mirrored or cancelled, as one deems appropriate. The only drawback would be the resulting increase in phase lag and group delay.
RHP zero => Mirror to LHP pole (flat magnitude, double phase)
LHP zero => Cancel using LHP pole (flat magnitude, zero phase)
Nevertheless, Mr. Berchin's z-plane method's still there. The MATLAB code is also quite old anyway, having been written in '86, with revisions in '88 and '92.
Last edited:
That's the mechanism to invoke iteration but I think the rationale is as I said.
Yes, that's my kind of approach, simple and obvious has an appeal!
However, the room correction systems don't usually try this, AFAIK.
Some of them just leave RHP zeroes alone I believe.
If I understand correctly, that "Invfreqs" tries to find a pole zero pattern that is close to the specified response but with a stable inverse, then that is new to me.
Thanks for your contribution of the documentation for this.
Best wishes
David