[speex-dev] Fixed-point cos/acos

[Mike Dunn]

Jean-Marc Valin wrote:

>Hi,
>
>Before I try to code this myself, I'd like to know if anyone has a
>fixed-point routine to compute the cos and acos functions. All I need is
>around 3-digit accuracy. 
>
>Thanks.
>
>	Jean-Marc
>
>  
>
I think I can help with cos.  icosine() and isine() below take input in 
"bnary angular measure", or "pirad" format, where the full range 
represents -pi to +pi (or equivalently, 0 to 2*pi for unsigned shorts). 
 The input is translated into the +/- pi/4 range by isine(), and using 
the trig identities, the appropriate _sin or _cos function is called. 
 These functions perform the series expansion with three terms, which 
should provide the accuracy you need given their constrained inputs.

BTW, this is based on the material in Jack Crenshaw's excellent book 
"Math Toolkit for Real-Time Programming".

HTH,
Mike Dunn

hort _sin(short y)
{
  /* Calculate sin using series expansion with 3 terms.
     Input of this function is limited to +/- pi/4 */

  /* series coefficients, tweaked with Chebyshev */
  const short s1 = 0x6487;
  const short s3 = 0x2951;
  const short s5 = 0x4f6;

  long z, prod, sum;

  z = ((long)y * y) >> 12;
  prod = (z * s5) >> 16;
  sum = s3 - prod;
  prod = (z * sum) >> 16;
  sum = s1 - prod;

  return (short)((y * sum) >> 13);
}

hort _cos(short y)
{
  const short c0 = 0x7fff;
  const short c2 = 0x4eea;
  const short c4 = 0x0fc4;

  long z, prod, sum;
  z = ((long)y * y) >> 12;
  prod = (z * c4) >> 16;
  sum = c2 - prod;

  prod = (z * sum) >> 15;
  return (short)(c0 - prod);
}

<p>short isine(short x)
{
  /* Input of this function is in binary angular measure format (Q15
     with the full -1 to 1 range representing -pi to pi). */

  /* Determine the quadrant the input resides in */
  unsigned short n = ((unsigned short)x + 0x2000) >> 14;
  x -= n * 0x4000;     /* translate input down to +/- pi/4 */

  /* call the appropriate function given the quadrant of the input */
  switch(n){
  case 0:
    return _sin(x);
  case 1:
    return _cos(x);
  case 2:
    return - _sin(x);
  case 3:
    return  - _cos(x);
  }
}

<p>short icosine(short x){
  return isine(x + 0x4000);
}

<p><p><p><p><p><p>--- >8 ----
List archives:  http://www.xiph.org/archives/
Ogg project homepage: http://www.xiph.org/ogg/
To unsubscribe from this list, send a message to 'speex-dev-request at xiph.org'
containing only the word 'unsubscribe' in the body.  No subject is needed.
Unsubscribe messages sent to the list will be ignored/filtered.