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