4 ms·

I recently implemented the exponential function for a sound synthesizer toolkit I’m working on. The method I used was not mentioned in the article, so I’ll expl

I recently implemented the exponential function for a sound synthesizer toolkit I’m working on. The method I used was not mentioned in the article, so I’ll explain it here. I used the Remez algorithm, modified to minimize equivalent input error. This algorithm lets you find a polynomial with the smallest maximum error. You start with a set of X coordinates, and create a polynomial which oscillates up and down around the target function, so p(x0) = f(x0) + E, p(x1) = f(x1) - E, p(x2) = f(x2) + E, etc. You then iterate by replacing the set of x values with the x values which have maximum error. This will converge to a polynomial with smallest maximum error. From my experiments, you can use a fairly low order approximation for musical applications. Used to calculate frequencies of musical notes, 2nd order is within 3.0 cents, and 3rd order is within 0.13 cents.

exmadscientist 6y ago
Remez is the most common implementation of Minimax approximation algorithms, so it's probably coming when the author makes it to section 6.
- fractionalhare 6y ago
  Yes, this is the plan :)
  - deleted 6y ago
    [deleted]
dvt 6y ago
This is way cool! Mad respect to you audio engineers out there. I tried implementing FFT last year in one of my open-source projects[1] (key word: tried). [1] https://github.com/dvx/lofi/commit/285bf80ff6d7f0784c14270de487f96e076ffca6 https://github.com/dvx/lofi/commit/285bf80ff6d7f0784c14270de...
- nayuki 6y ago
  Glancing at your code diff, I think you're missing the bit-reversed permutation. Either way, feel free to steal from: https://www.nayuki.io/page/free-small-fft-in-multiple-languages https://www.nayuki.io/page/free-small-fft-in-multiple-langua...
jbay808 6y ago
I wouldn't recommend doing this directly -- there's no good polynomial that can represent the exponential function well over a wide range. Instead, it's better to exploit the definition of IEEE754 as an exponent and a mantissa. You calculate the exponent directly (because e^x = 2^(x/ln2)), and use the Remez algorithm to find a polynomial that fits just the mantissa.
- deleted 6y ago
  [deleted]
- dietrichepp 6y ago
  Honestly, reading this comment made me a bit unhappy. Don't mansplain to me the things that I've done. The Remez algorithm requires that you specify the domain. If you tried to use it to approximate the exponential function over its entire domain the algorithm would diverge. You have to explicitly choose a range [x1,x2] for the Remez algorithm. It minimizes the maximum error over that range. I'm ordinarily happy to share the work I'm doing but this comment reminds me of the reasons I don't share much with HN.
  - 6nf 6y ago
    > Don't mansplain to me I honestly don't think this has anything to do with your gender or the gender of the person you're replying to.
  - jbay808 6y ago
    I'm writing this because I've done it too, for a performance-critical embedded application, and I think it's an insight worth sharing, for anybody else who comes across this and is trying to do something similar. Not to put you down or one-up you or anything; I'm very sorry to have come across that way. Can I try again? For a domain like [0,1], the Remez algorithm will do a great job -- because the exponent is barely changing, just the mantissa. For a range like [-10, 10], you won't find any polynomial of reasonable order that has acceptable error. But, there's a really neat trick that does let the Remez algorithm work over the entire domain of floats, and get excellent performance with just a fourth or fifth order polynomial, or even third order if you're pressed for time, and minimal extra computation. That is, with just a little simplification, you multiply x by (1/ln(2)), floor it, and stuff the resulting integer into the exponent bits. Then, you take the remainder (what is left over after taking the floor), and make that your input to the polynomial, and stuff the result into the mantissa bits. There's a little more to it; signs and NaNs need to be handled correctly, subnormal numbers can be treated as a special case with a different Remez polynomial, and a few other corner cases and exact values might be important (e^0, e^1, e^-1). But the result works wonders and it's basically just using Remez but transforming the input to a domain that behaves more like a polynomial, leveraging the magic of IEEE-754.
nitrogen 6y ago
I used something similar to implement an approximate square root function in a graphical DSP platform that only has addition and multiplication. A 4th-order polynomial gives satisfactory accuracy from -20dBFS to 0 dBFS, which was the region I cared about. It's hard to do better near the origin (-inf dB) though, because the slope of the square root approaches infinity.

exmadscientist 6y ago

Remez is the most common implementation of Minimax approximation algorithms, so it's probably coming when the author makes it to section 6.

fractionalhare 6y ago
Yes, this is the plan :)
- deleted 6y ago
  [deleted]

dvt 6y ago

This is way cool! Mad respect to you audio engineers out there. I tried implementing FFT last year in one of my open-source projects[1] (key word: tried). [1] https://github.com/dvx/lofi/commit/285bf80ff6d7f0784c14270de487f96e076ffca6 https://github.com/dvx/lofi/commit/285bf80ff6d7f0784c14270de...

nayuki 6y ago
Glancing at your code diff, I think you're missing the bit-reversed permutation. Either way, feel free to steal from: https://www.nayuki.io/page/free-small-fft-in-multiple-languages https://www.nayuki.io/page/free-small-fft-in-multiple-langua...

jbay808 6y ago

I wouldn't recommend doing this directly -- there's no good polynomial that can represent the exponential function well over a wide range. Instead, it's better to exploit the definition of IEEE754 as an exponent and a mantissa. You calculate the exponent directly (because e^x = 2^(x/ln2)), and use the Remez algorithm to find a polynomial that fits just the mantissa.

deleted 6y ago
[deleted]
dietrichepp 6y ago
Honestly, reading this comment made me a bit unhappy. Don't mansplain to me the things that I've done. The Remez algorithm requires that you specify the domain. If you tried to use it to approximate the exponential function over its entire domain the algorithm would diverge. You have to explicitly choose a range [x1,x2] for the Remez algorithm. It minimizes the maximum error over that range. I'm ordinarily happy to share the work I'm doing but this comment reminds me of the reasons I don't share much with HN.
- 6nf 6y ago
  > Don't mansplain to me I honestly don't think this has anything to do with your gender or the gender of the person you're replying to.
- jbay808 6y ago
  I'm writing this because I've done it too, for a performance-critical embedded application, and I think it's an insight worth sharing, for anybody else who comes across this and is trying to do something similar. Not to put you down or one-up you or anything; I'm very sorry to have come across that way. Can I try again? For a domain like [0,1], the Remez algorithm will do a great job -- because the exponent is barely changing, just the mantissa. For a range like [-10, 10], you won't find any polynomial of reasonable order that has acceptable error. But, there's a really neat trick that does let the Remez algorithm work over the entire domain of floats, and get excellent performance with just a fourth or fifth order polynomial, or even third order if you're pressed for time, and minimal extra computation. That is, with just a little simplification, you multiply x by (1/ln(2)), floor it, and stuff the resulting integer into the exponent bits. Then, you take the remainder (what is left over after taking the floor), and make that your input to the polynomial, and stuff the result into the mantissa bits. There's a little more to it; signs and NaNs need to be handled correctly, subnormal numbers can be treated as a special case with a different Remez polynomial, and a few other corner cases and exact values might be important (e^0, e^1, e^-1). But the result works wonders and it's basically just using Remez but transforming the input to a domain that behaves more like a polynomial, leveraging the magic of IEEE-754.

nitrogen 6y ago

I used something similar to implement an approximate square root function in a graphical DSP platform that only has addition and multiplication. A 4th-order polynomial gives satisfactory accuracy from -20dBFS to 0 dBFS, which was the region I cared about. It's hard to do better near the origin (-inf dB) though, because the slope of the square root approaches infinity.