Here is an example of a real-life problem I needed to solve, using NML. But first, why do I refer to NML as the Ultimate Lambda for Scientific Programming?
I don’t actually believe that claim. The Ultimate Ultimate Language for any sort of programming, scientific or otherwise, is undoubtedly Lisp. And I use Lisp almost as often as I use NML. Lisp is the ultimate power weapon of the universe and its macrology enables the production of application specific languages in a way that the fixed-syntax language NML could not hope to duplicate.
But Lisp is also very cumbersome at times, and NML is specifically optimized for application to vectorized math problems. NML is very very terse compared to Lisp for numerical problems, and so I almost always turn to NML first for quick prototyping.
So here is the real-world problem I was facing... I needed a simple rational approximation to a highly nonlinear function related to nonlinear gain compression in an audio system.
The nonlinear gain can only be computed as the root of a transcendental equation. That requires an iterative solution, entailing some unknown number of iterations. Hence the time required to compute the nonlinear gain would depend on the instantaneous input loudness and some lucky guess as to the starting estimate for the answer.
The gain compression function needed to be as accurate as possible over the widest typical range of audio levels. And furthermore, it needed to maintain that accuracy under a nonlinear weighting that exponentially increased toward the lower audio levels.
So I wrote a routine for producing MiniMax rational approximations with arbitrary weighting functions. I needed a simple approximation so that the nonlinear compression gain could be computed rapidly in a real-time DSP environment.
Simple polynomial approximations could be used, but to achieve the desired accuracy, that would have meant very high order polynomials under some circumstances. So I settled for a rational approximation which incurs a divide operation (which my DSP - an Intel Pentium IV has) and lower net order in the numerator and denominator polynomials for a given level of accuracy.
I settled on an order (2,2) rational approximation as most suitable for my range of control parameters. And here is an example of the minimax error for such a rational approximation at one particular control parameter level:
So under the weighting shown above this rational approximation promises a maximum absolute error of less than 0.15 dB over the entire input loudness domain. Our x-axis in the error plot runs from -1 to +1 as a result of our use of a range-reduction equation that maps the input range of [20, 100] phons into the approximator range of [-1, +1].
What does this approximation actually look like? Here it is:
Gdb = (5.06176561961-9.67731394114x+5.5277831969x^2) /
(1+0.624093143413x+0.281860954703x^2)
The code starts out by assuming the nodes in the error curve occur at the nodes of a Tchebyshev polynomial of degree (Np + Nq + 1) where Np is the order of the numerator polynomial, and Nq is the order of the denominator polynomial. Of course, under the highly nonlinear weighting we used, that won’t be correct, but it ought to get us in the neighborhood of where we can start to iterate toward the correct nodes.
However, under some circumstances, those starting nodes will not be near enough to the correct ones to allow successful convergence to the actual nodes. In that case, we raise an exception in the code, and trap that exception in a small routine that will randomly jiggle the locations of the nodes and then retry the fitting routine. If we fail to converge after some number of restarts, then we simply give up and propagate the exception to the caller of the Remez fitting routine.
Included in the code are routines for locating the extrema of the error curve using some heuristics along with narrowed interval bisection. We found Newton’s method or quadratic interpolation for locating the extrema to be unreliable under some circumstances. And since we don’t have analytic functions for the derivative of the error curve, we wing it using interval bisection. The result is still blazing fast.
We compared the speed of our approximation routine, which takes time out to visually plot the progress of the fitting, to a Mathematica routine for solving the same problem using their built-in rational minimax fitting function. Our solution in NML was over 100 times faster.
And since we have the source code for the fitting routines, we can modify them as needed, insert graphing diagnostics and progress indications, and do anything else that we want. There is no limit to our ability to extend this solution to other areas where some modification would be needed. That isn’t true of the built-in Mathematica routine.
Saturday, February 3, 2007
