6
$\begingroup$

I'm looking for books or resources that cover the following in detail:

  • implementing mathematical functions (e.g., logarithm, exponential, sine, cosine, inverse) in fixed point arithmetic for DSP purposes.

  • techniques like using lookup tables, Taylor series, etc.

I am fairly familiar with C programming and am more interested in the algorithms on how to go about implementing various mathematical functions in an efficient manner.

$\endgroup$
  • 1
    $\begingroup$ This is just one trick, but it's very useful. It's about computing the atan2 function, i.e. computing the argument of a complex number from its real and imaginary parts. $\endgroup$ – Matt L. Feb 10 '15 at 11:31
  • 1
    $\begingroup$ i can give you some optimized finite-order power series that i developed more than a decade ago. other than the initial client, i haven't gotten any money for it since so i think i may as well make it public domain. originally it was developed for floating point in a context where the client couldn't include the stdlib in the build. and the evaluation is not bit perfect. i.e. there is error, but very small and optimized for the particular function being evaluated. lemme find that file, i'll pull the coefficients out and post the series's. $\endgroup$ – robert bristow-johnson Feb 11 '15 at 0:43
  • $\begingroup$ @robertbristow-johnson i look forward to making use of the series and seeing how it goes! thanks! $\endgroup$ – RuD Feb 11 '15 at 7:18
  • $\begingroup$ Ruchir, i posted the series's below. there are common-sense things you need to do to extend the range. like $$ 2^x = 2^k \times 2^{x-k} $$ if $ k \le x \le k+1 $. similar thing for $\log_2()$ and the periodic sinusoids. for exp and log, this will require arithmetic bit shifting of what comes out (for exp) or what goes in (for log). i think you can figger that out. and, of course, for exp and log of different bases (like $e$), you just scale with the appropriate constant what goes in to $2^x$ and what comes outa $\log_2(x)$. $\endgroup$ – robert bristow-johnson Feb 12 '15 at 16:26
8
$\begingroup$

the general polynomial form is:

$$\begin{align} f(u) &= \sum\limits_{n=0}^{N} \ a_n \ u^n \\ \\ &= a_{\small{0}} + \Bigg(a_{\small{1}} + \bigg(a_{\small{2}} + \Big(a_{\small{3}} + \,... \big(a_{\small{N-2}} + (a_{\small{N-1}} + a_{\small{N}} \,u \,)u \, \big)u \ ...\Big)u \, \bigg)u \, \Bigg)u\\ \end{align}$$

the latter form is using Horner's method, which is highly recommended, especially if you're doing this in single-precision floating point.

then for a few specific functions:

square root:

$$ \begin{align} f(x-1) & \approx \sqrt{x} \quad \quad 1 \le x \le 2 \quad \quad N=4\\ a_0 & = 1.0 \\ a_1 & = 0.49959804148061 \\ a_2 & = -0.12047308243453 \\ a_3 & = 0.04585425015501 \\ a_4 & = -0.01076564682800 \\ \end{align} $$

if $2 \le x \le 4$, use the above to evaluate $\sqrt{\tfrac{x}{2}}$ and multiply that result with $\sqrt{2}$ to get $\sqrt{x}$. as with $\log_2(x)$, apply power of $2$ scaling to scale the argument to the necessary range.

base 2 logarithm:

$$ \begin{align} x\cdot f(x-1) & \approx \log_2(x) \quad \quad 1 \le x \le 2 \quad \quad N=5\\ a_0 & = 1.44254494359510 \\ a_1 & = -0.7181452567504 \\ a_2 & = 0.45754919692582 \\ a_3 & = -0.27790534462866 \\ a_4 & = 0.121797910687826 \\ a_5 & = -0.02584144982967 \\ \end{align} $$

base 2 exponential:

$$ \begin{align} f(x) & \approx 2^x \quad \quad 0 \le x \le 1 \quad \quad N=4\\ a_0 & = 1.0 \\ a_1 & = 0.69303212081966 \\ a_2 & = 0.24137976293709 \\ a_3 & = 0.05203236900844 \\ a_4 & = 0.01355574723481 \\ \end{align} $$

sine:

$$ \begin{align} x\cdot f(x^2) & \approx \sin\left(\tfrac{\pi}{2} x \right) \quad \quad -1 \le x \le 1 \quad \quad N=4 \\ a_0 & = 1.57079632679490 \\ a_1 & = -0.64596406188166 \\ a_2 & = 0.07969158490912 \\ a_3 & = -0.00467687997706 \\ a_4 & = 0.00015303015470 \\ \end{align} $$

cosine (use sine):

$$ \cos(\pi x) = 1 \, - \, 2 \, \sin^2 \left(\tfrac{\pi}{2} x \right) $$

tangent:

$$ \tan(x) = \frac{\sin(x)}{\cos(x)} $$

inverse tangent:

$$ \begin{align} \frac{x}{f(x^2)} & \approx \arctan(x) \quad \quad -1 \le x \le 1 \quad \quad N=4 \\ a_0 & = 1.0 \\ a_1 & = 0.33288950512027 \\ a_2 & = -0.08467922817644 \\ a_3 & = 0.03252232640125 \\ a_4 & = -0.00749305860992 \\ \end{align} $$

$$ \arctan(x) = \tfrac{\pi}{2} - \arctan\left( \tfrac{1}{x} \right) \quad \quad 1 \le x $$

$$ \arctan(x) = -\tfrac{\pi}{2} - \arctan\left( \tfrac{1}{x} \right) \quad \quad x \le -1 $$

inverse sine:

$$ \arcsin(x) = \arctan\left( \frac{x}{\sqrt{1-x^2}} \right)$$

inverse cosine:

$$\begin{align} \arccos(x) &= \frac{\pi}{2} - \arcsin(x) \\ &= \frac{\pi}{2} - \arctan\left( \frac{x}{\sqrt{1-x^2}} \right)\\ \end{align}$$

$\endgroup$
  • $\begingroup$ That seems pretty useful! Thank you a lot for sharing that. But still first reference entry in man sox is best ;) $\endgroup$ – jojek Feb 12 '15 at 15:52
  • $\begingroup$ dunno sox. what does the manual say about it? $\endgroup$ – robert bristow-johnson Feb 12 '15 at 15:56
  • 2
    $\begingroup$ Simply [1] R. Bristow-Johnson, Cookbook formulae for audio EQ biquad filter coefficients, http://musicdsp.org/files/Audio-EQ-Cookbook.txt :) $\endgroup$ – jojek Feb 12 '15 at 15:58
  • $\begingroup$ BTW, the series's minimizes the maximum weighted error. error is weighted in such a manner that makes sense for the function. uniform maximum error for $\log_2()$. proportional maximum error for $\sqrt{x}$ and $2^x$. something similar to proportional for $\sin()$ having to do with computing coefficients for resonant filters so that the max error of $\log(f_0)$ is minimized. i can't remember what criterion i used for $\arctan()$. and for some reason i can't find my file that tells me what the maximum errors are, given the range of $x$. someone with MATLAB can find out. $\endgroup$ – robert bristow-johnson Feb 12 '15 at 16:04
  • $\begingroup$ the other thing i tried to do was to have zero error at the end-points (so the remez exchange alg was modified a little), and i chose the order of the series (even or odd) so that when the functions are spliced with copies of themselves to extend the range, the slope of the error is nearly continuous at the splice. if you plot the difference of the approximation to the real thing, you'll see what i mean. $\endgroup$ – robert bristow-johnson Feb 12 '15 at 16:10
2
$\begingroup$

Although not specific to fixed point, I would highly recommend the "Math Toolkit for Real-Time Programming" book by Jack Crenshaw. It comes with a CD with the source code.

$\endgroup$
2
$\begingroup$

TI has IQMath libraries for all of their fixed point microcontrollers. I have found them to be a goldmine of fixed-point math and DSP functions not necessarily limited to TI chips.

MSP430 C28X

$\endgroup$
  • $\begingroup$ I'm more interested in the algorithms rather than just the implementation of the functions $\endgroup$ – RuD Feb 16 '15 at 14:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.