Books/resources for implementing various mathematical functions in fixed point arithmetic for DSP purposes

Question

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.

Answer 1

C code that does this here

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}$$

Answer 2

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.

Answer 3

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

Answer 4

Chebyshev approximation can help compute polynomial coefficients that are close to optimum for approximating a function over a finite range. You run the approximation routine on a PC to achieve a particular set of polynomial coefficients, which you can then apply on whatever platform you like (e.g. embedded / DSP) The fine print is more or less as follows:

• This only works for functions of one variable; if you have some function $z=f(x,y)$ then Chebyshev approximation isn't going to help you.
• The function you're approximating should be "polynomial-like". Corners, sharp bends, and lots of wiggles will require higher-order polynomials to achieve a given level of accuracy.
• Keeping the polynomial order low is important -- above 5th or so, you may start to see numerical errors.
• Chebyshev approximation uses the Chebyshev polynomials evaluated over the domain [-1,1], so if the input range for your function is significantly different, you may need to scale/offset the input appropriately before determining the coefficients and before applying them. For example, if I care about a function over the input range $x \in [0,20]$ then I might define $u = (x\times 0.1) - 1$ so that $u$ ranges from -1 to +1 and I can evaluate a polynomial in $u$ to compute the required result. Without this scaling, you can run into numerical errors more easily -- or stated another way, for the same precision you may need higher bit lengths to compute intermediate values, and this is usually undesirable.