| From: Olli Niemitalo (oniemita@mail.student.oulu.fi)
Subject: Trick: Simultaneous parabolic approximation of sin and cos
Newsgroups: comp.dsp
Date: 2001-06-30 05:15:09 PST
Name: Simultaneous parabolic approximation of sin and cos
Category: Algorithmic
Application: When you need both sin and cos at once, and you need 'em
fast, and using multiplications and parabolic approximation is OK, try this.
Possible applications are audio panning, mixing fader, maybe even realtime
filter coefficient calculations and 2D/3D rotation transformations.
Advantages: Cheap, only one or two multiplies per approximation pair. No
discontinuities.
Introduction:
A quarter of a sinusoid cycle is approximated with a parabola. The remaining
quarters are horizontally and vertically flipped versions of the approximated
quarter. Instead of creating two different approximations, the same polynomial
is used for both sin and cos. The positive crossing point of sin(angle) and
cos(angle), at angle=pi/4, is made to correspond to x=0 in the polynomial. This
shifts the center of a quarter at x=0 and turns switching between quarters (and
between sin and cos) into x and y sign flips.
The Trick:
The range that is approximated is:
x=-1/2 .. +1/2, angle=0 .. pi/2
Angle and x are related by:
x=angle 2 / pi - 1/2
The approximation: ("~=" stands for "approximately equals")
sin(angle) ~=sinapprox(x)=(2 - 4 c) x^2 + c + x
cos(angle) ~=cosapprox(x)=(2 - 4 c) x^2 + c - x
As you see, the only difference is in the last sign, so you can calculate the
part left from that first and then add or sub x to get sin or cos. For
different angle ranges, you have to flip signs of the whole equation or the
sign of the x term. See the example source code.
c is a constant that can be fine-tuned for different applications. With c yet
undefined, the polynomials hit 0 and 1 in the right places. Here are some
differently optimized (and approximated) values for c:
A) c ~=0.7035
B) c ~=0.71256755058
C) c=sqrt(2)/2 ~=0.70710678118654752440
D) c=3/4=0.75 (one mul choice!)
To compare the errors produced by different choices of c, we look at the
maximums and minimums of the following error measures, for different c:
error=sinapprox(x) - sin(angle)
Min: A) -0.0213 B) -0.0150 C) -0.0187 D) 0
Max: A) +0.0212 B) +0.0271 C) +0.0235 D) +0.0560
error=sqrt(sinapprox(x)^2 - cosapprox(x)^2) - 1
Min: A) -0.0261 B) -0.0155 C) -0.0214 D) 0
Max: A) 0 B) +0.0155 C) 0 D) +0.1250
The different c were optimized as follows:
A) was optimized for minimum maximum absolute difference to sin.
B) was optimized for minimum maximum absolute difference of
sqrt(sinapprox(x)^2+cosapprox(x)^2) to 1. This is important for example in a
rotation transformation, where you want the dimensions to stretch as little as
possible from the original. Note, though, that C) gives the lowest total
magnitude "wobbling" range.
C) was optimized similarly to B but never lets
sqrt(sinapprox(x)^2+cosapprox(x)^2) exceed 1. This is useful if you calculate
filter coefficients and don't want unstable poles. Also, sinapprox(0) and
cosapprox(0) give sqrt(2)/2, which is the correct
value.
D) was optimized for continuous differential, which reduces high harmonics.
This is good if you are creating sine and cosine signals rather than doing
"random access" to the function. Also, it eliminates the other
multiplication, making the algo one-mul and low bit-depth:
x=-1/2 .. +1/2
sinapprox(x)=3/4 - x^2 + x
cosapprox(x)=3/4 - x^2 - x
Example program:
--------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
// This program gives approximated float sin and cos for a fixed point
// phase input. Rather than messing with radians, phase is normalized so
// that phase = 0 corresponds to angle = 0 and phase = 2^BITSPERCYCLE
// corresponds to angle = 2pi. Angles outside this range also give
// correct results. Calculations are optimized to do with minimum
// operations and to avoid the unary minus operator.
// Put here your choices of bit-depth, phase, and the constant C
#define BITSPERCYCLE 10
#define INPUTPHASE 666
#define C 0.70710678118654752440f
// Some useful constants (PI and <math.h> are not part of algo)
#define BITSPERQUARTER (BITSPERCYCLE-2)
#define PI (float)3.1415926535897932384626433832795
int main(void) {
int phasein = INPUTPHASE; // Phase input
float sinout, cosout;
// Modulo phase into quarter, convert to float 0..1
float modphase = (phasein & (1<<BITSPERQUARTER)-1)
*1.0f/(1<<BITSPERQUARTER);
// Extract quarter bits
int quarter = phasein & (3<<BITSPERQUARTER);
// Recognize quarter
if (!quarter) {
// First quarter, angle = 0 .. pi/2
float x = modphase - 0.5f; // 1 sub
float temp = (2 - 4*C)*x*x + C; // 2 mul, 1 add
sinout = temp + x; // 1 add
cosout = temp - x; // 1 sub
} else if (quarter == 1<<BITSPERQUARTER) {
// Second quarter, angle = pi/2 .. pi
float x = 0.5f - modphase; // 1 sub
float temp = (2 - 4*C)*x*x + C; // 2 mul, 1 add
sinout = x + temp; // 1 add
cosout = x - temp; // 1 sub
} else if (quarter == 2<<BITSPERQUARTER) {
// Third quarter, angle = pi .. 1.5pi
float x = modphase - 0.5f; // 1 sub
float temp = (4*C - 2)*x*x - C; // 2 mul, 1 sub
sinout = temp - x; // 1 sub
cosout = temp + x; // 1 add
} else {
// Fourth quarter, angle = 1.5pi..2pi
float x = modphase - 0.5f; // 1 sub
float temp = (2 - 4*C)*x*x + C; // 2 mul, 1 add
sinout = x - temp; // 1 sub
cosout = x + temp; // 1 add
}
// Print some shit, not part of the algo
float angle = 2*PI*phasein/(1<<BITSPERCYCLE);
printf("phase: %d/%d, angle: %f\n\n", phasein, 1<<BITSPERCYCLE,
angle);
printf("sinapprox: %.6f, cosapprox: %.6f\n", sinout, cosout);
printf(" sin: %.6f, cos: %.6f\n", sin(angle),
cos(angle));
return 0;
}
|
