Taylor series equation - angle value setup

susan-parker

Member
I have tried this out with a certain new RA4M1 board, but found that the "angle" setup needed to be:

angle = (int32_t)(0x7FFFFFFFu - ph);

... otherwise at the very top point of the sine these was a glitch (which I noticed as I was using the DAC as a monitor).

This was using an increment derived from:

uint32_t sine_fraction_delta = 0x1 << 24;

to step through the 32 bit uint value.

sine_phase_accumulate += sine_fraction_delta ;

Original code details from here:

https://www.pjrc.com/high-precision-sine-wave-synthesis-using-taylor-series/

// High accuracy 11th order Taylor Series Approximation
// input is 0 to 0xFFFFFFFF, representing 0 to 360 degree phase
// output is 32 bit signed integer, top 25 bits should be very good
static int32_t taylor(uint32_t ph)
{
int32_t angle, sum, p1, p2, p3, p5, p7, p9, p11;

if (ph >= 0xC0000000 || ph < 0x40000000) {
angle = (int32_t)ph; // valid from -90 to +90 degrees
} else {
angle = (int32_t)(0x80000000u - ph);
}
p1 = multiply_32x32_rshift32_rounded(angle << 1, 1686629713);
p2 = multiply_32x32_rshift32_rounded(p1, p1) << 3;
p3 = multiply_32x32_rshift32_rounded(p2, p1) << 3;
sum = multiply_subtract_32x32_rshift32_rounded(p1 << 1, p3, 1431655765);
p5 = multiply_32x32_rshift32_rounded(p3, p2) << 1;
sum = multiply_accumulate_32x32_rshift32_rounded(sum, p5, 286331153);
p7 = multiply_32x32_rshift32_rounded(p5, p2);
sum = multiply_subtract_32x32_rshift32_rounded(sum, p7, 54539267);
p9 = multiply_32x32_rshift32_rounded(p7, p2);
sum = multiply_accumulate_32x32_rshift32_rounded(sum, p9, 6059919);
p11 = multiply_32x32_rshift32_rounded(p9, p2);
sum = multiply_subtract_32x32_rshift32_rounded(sum, p11, 440721);
return sum <<= 1;
}

And it is wonderfully fast for such an otherwise "slow" processor.
Thanks.

I agree: taylor(1073741824) returns = -2147483522. When lines 84-88 of https://github.com/PaulStoffregen/Audio/synth_sine.cpp are changed to
Code:
``````	if (ph >= 0xC0000000 || ph < 0x40000000) {                            // ph:  0.32
angle = (int32_t)ph; // valid from -90 to +90 degrees
} else {
angle = (int32_t)(0x7FFFFFFFu - ph);                        // angle: 2.30
}``````
it will return 2147483524 instead; the exact correct value would be 2147483647. This limits the absolute error in taylor() to about -129 .. +135 (using an exhaustive comparison for all possible phases, compared to double-precision values).

Many thanks for the confirm...

I agree: taylor(1073741824) returns = -2147483522.

Many thanks for the confirm.

I have just tidied up my test-code and checked it runs on a Teensy 4:

Code:
``````/*  Teensy 4 SINE wave output glitch diagnosis
*
*  From RA4M1 DAC code
*
*  Susan Parker - 21st August 2023.
*    Undate with Paul Stoffregen's 11th order Taylor Series Approximation
*    https://github.com/PaulStoffregen/Audio/blob/master/utility/dspinst.h
*
*  Susan Parker - 22nd August 2023.
*    Debug Taylor Series Approximation glitch at +90.000° = 0x40000000
*/

// From "Quick loop DAC output test - with sine/cosine calls"

#define SINE_GLITCH_DIAG  // Print out values around problem state
#define SINE_GLITCH_LOOP  // Use 0x80000000 - otherwise 0x7FFFFFFF

#include <Arduino.h>
#include "synth_sine.h"
#include "utility/dspinst.h"

#ifdef SINE_GLITCH_DIAG
int32_t angle;
#endif

void setup()
{
int32_t result_taylor         = 0;
uint32_t sine_phase_accumulate = 0;

Serial.begin(115200);      // The interrupts for the USB serial are already in place before setup() starts
while (!Serial){};         // Note: USB serial cannot be used for serial comms when running fast IRQs - diagnostics only.

#ifdef SINE_GLITCH_DIAG
Serial.println("Fault at 0x40000000");
for(sine_phase_accumulate = 0x40000000-8; sine_phase_accumulate < 0x40000000+8; sine_phase_accumulate++)
{
result_taylor = taylor_line(sine_phase_accumulate);
Serial.print(sine_phase_accumulate, HEX);
Serial.print(" - ");
Serial.print(result_taylor, HEX);
Serial.print(" - ");
Serial.println(angle, HEX);
}

Serial.println("No Fault at 0xC0000000");
for(sine_phase_accumulate = 0xC0000000-8; sine_phase_accumulate < 0xC0000000+8; sine_phase_accumulate++)
{
result_taylor = taylor_line(sine_phase_accumulate);
Serial.print(sine_phase_accumulate, HEX);
Serial.print(" - ");
Serial.print(result_taylor, HEX);
Serial.print(" - ");
Serial.println(angle, HEX);
}
#endif
}

void loop(void)                        //
{                                    //
delay(1);
}

// High accuracy 11th order Taylor Series Approximation
// input is 0 to 0xFFFFFFFF, representing 0 to 360 degree phase
// output is 32 bit signed integer, top 25 bits should be very good

static int32_t taylor_line(uint32_t ph)  // Inline version - this code is to go inside an interrupt call
{
#ifdef SINE_GLITCH_DIAG
int32_t sum, p1, p2, p3, p5, p7, p9, p11, term;
#else
int32_t angle, sum, p1, p2, p3, p5, p7, p9, p11, term;
#endif

if (ph >= 0xC0000000u || ph < 0x40000000u)
{
angle = (int32_t)ph; // valid from -90 to +90 degrees
}
else
{
#ifdef SINE_GLITCH_LOOP
angle = (int32_t)(0x80000000u - ph);
#else
angle = (int32_t)(0x7FFFFFFFu - ph);
#endif
}
term = angle << 1;
asm volatile("smmulr %0, %1, %2" : "=r" (p1) : "r" (term), "r" (1686629713));
asm volatile("smmulr %0, %1, %2" : "=r" (term) : "r" (p1), "r" (p1));
p2 = term << 3;
asm volatile("smmulr %0, %1, %2" : "=r" (term) : "r" (p2), "r" (p1));
p3 = term << 3;
term = p1 << 1;
asm volatile("smmlsr %0, %2, %3, %1" : "=r" (sum) : "r" (term), "r" (p3), "r" (1431655765));
asm volatile("smmulr %0, %1, %2" : "=r" (term) : "r" (p3), "r" (p2));
p5 = term << 1;
asm volatile("smmlar %0, %2, %3, %1" : "=r" (sum) : "r" (sum), "r" (p5), "r" (286331153));
asm volatile("smmulr %0, %1, %2" : "=r" (p7) : "r" (p5), "r" (p2));
asm volatile("smmlsr %0, %2, %3, %1" : "=r" (sum) : "r" (sum), "r" (p7), "r" (54539267));
asm volatile("smmulr %0, %1, %2" : "=r" (p9) : "r" (p7), "r" (p2));
asm volatile("smmlar %0, %2, %3, %1" : "=r" (sum) : "r" (sum), "r" (p9), "r" (6059919));
asm volatile("smmulr %0, %1, %2" : "=r" (p11) : "r" (p9), "r" (p2));
asm volatile("smmlsr %0, %2, %3, %1" : "=r" (sum) : "r" (sum), "r" (p11), "r" (440721));
return sum <<= 1;
}``````

I am super impressed with the speed and accuracy of this maths, which is way beyond my capacity to do from scratch. Just sorting out the glitch took me several hours!

Last edited:
This is what the glitch shows as in the test-code output:

Fault at 0x40000000
3FFFFFF8 - 7FFFFF88 - 3FFFFFF8
3FFFFFF9 - 7FFFFF84 - 3FFFFFF9
3FFFFFFA - 7FFFFF88 - 3FFFFFFA
3FFFFFFB - 7FFFFF88 - 3FFFFFFB
3FFFFFFC - 7FFFFF88 - 3FFFFFFC
3FFFFFFD - 7FFFFF8A - 3FFFFFFD
3FFFFFFE - 7FFFFF88 - 3FFFFFFE
3FFFFFFF - 7FFFFF84 - 3FFFFFFF
40000000 - 8000007C - 40000000
40000001 - 7FFFFF84 - 3FFFFFFF
40000002 - 7FFFFF88 - 3FFFFFFE
40000003 - 7FFFFF8A - 3FFFFFFD
40000004 - 7FFFFF88 - 3FFFFFFC
40000005 - 7FFFFF88 - 3FFFFFFB
40000006 - 7FFFFF88 - 3FFFFFFA
40000007 - 7FFFFF84 - 3FFFFFF9
No Fault at 0xC0000000
BFFFFFF8 - 80000078 - C0000008
BFFFFFF9 - 8000007C - C0000007
BFFFFFFA - 80000078 - C0000006
BFFFFFFB - 80000078 - C0000005
BFFFFFFC - 80000078 - C0000004
BFFFFFFD - 80000076 - C0000003
BFFFFFFE - 80000078 - C0000002
BFFFFFFF - 8000007C - C0000001
C0000000 - 8000007C - C0000000
C0000001 - 8000007C - C0000001
C0000002 - 80000078 - C0000002
C0000003 - 80000076 - C0000003
C0000004 - 80000078 - C0000004
C0000005 - 80000078 - C0000005
C0000006 - 80000078 - C0000006
C0000007 - 8000007C - C0000007