Generate sine signal in C without using the standard function

Question

I want to generate a sine signal in C without using the standard function sin() in order to trigger sine shaped changes in the brightness of a LED. My basic idea was to use a lookup table with 40 points and interpolation.

Here's my first approach:

const int sine_table[40] = {0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767,  32364, 31163, 29196, 26509, 23170, 19260, 14876, 10125,
5125, 0, -5126, -10126,-14877, -19261, -23171, -26510, -29197, -31164, -32365,
-32768, -32365, -31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126};

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase)
{
    x1 = (int) phase % 41;
    x2 = x1 + 1;
    y = (sine_table[x2] - sine_table[x1])*((float) ((int) (40*0.001*i*100) % 4100)/100 - x1) + sine_table[x1];
    return y;
}

int main()
{
    while(1)
    {
    printf("%f      ", sin1(40*0.001*i)/32768);
    i = i + 1;
    }
}

Unfortunately, this function sometimes returns values far bigger than 1. Furthermore, the interpolation doesn't seem to be good (I used this to create sine shaped brightness changes of a LED, but these are very unsmoooth).

Does anybody have a better idea to implement a sine generator in C?

Answer 1

OP's main problem is in generating the index for the table look-up.

OP's code attempts to access outside array sine_table[40] leading to undefined behavior. Fix that at least.

const int sine_table[40] = {0, 5125, 10125, ...
    ...
    x1 = (int) phase % 41;                     // -40 <= x1 <= 40
    x2 = x1 + 1;                               // -39 <= x2 <= 41  
    y = (sine_table[x2] - sine_table[x1])*...  // bad code, consider x1 = 40 or x2 = 40,41

Suggested change

    x1 = (int) phase % 40;   // mod 40, not 41
    if (x1 < 0) x1 += 40;    // Handle negative values
    x2 = (x1 + 1) % 40;      // Handle wrap-around 
    y = (sine_table[x2] - sine_table[x1])*...  

There exist much better approaches, yet to focus on OP's method see below.

#include <math.h>
#include <stdio.h>

const int sine_table[40] = { 0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767, 32364, 31163, 29196, 26509, 23170, 19260, 14876, 10125,
5125, 0, -5126, -10126, -14877, -19261, -23171, -26510, -29197, -31164, -32365,
-32768, -32365, -31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126 };

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase) {
  x1 = (int) phase % 40;
  if (x1 < 0) x1 += 40;
  x2 = (x1 + 1) % 40;
  y = (sine_table[x2] - sine_table[x1])
      * ((float) ((int) (40 * 0.001 * i * 100) % 4100) / 100 - x1)
      + sine_table[x1];
  return y;
}

int main(void) {
  double pi = 3.1415926535897932384626433832795;
  for (int j = 0; j < 1000; j++) {
    float x = 40 * 0.001 * i;
    float radians = x * 2 * pi / 40;
    printf("%f %f %f\n", x, sin1(x) / 32768, sin(radians));
    i = i + 1;
  }
}

Output

         OP's     Reference sin()
0.000000 0.000000 0.000000
0.040000 0.006256 0.006283
0.080000 0.012512 0.012566
...
1.960000 0.301361 0.303035
2.000000 0.308990 0.309017
2.040000 0.314790 0.314987
...
39.880001 -0.020336 -0.018848
39.919998 -0.014079 -0.012567
39.959999 -0.006257 -0.006283

---

Better code would not pass the values i, x1, x2, y as global variables, but as function parameters or function variables. Perhaps that is an artifact of OP's debugging.

---

Does anybody have a better idea to implement a sine generator in C?

This is quite broad. Better as in speed, precision, code space, portability, or maintainability? sine() functions are easy to make. High-quality ones take more effort.

Although fuzzy, OP's use of a small look-up table is a good beginning - although I see it can be done without any floating point math. I recommend for OP to construct a tested and working solution and post it in Code Review for improvement ideas.

Answer 2

...a better idea to implement a sine generator in C?

Edit: Suggest first reading this article to gain an appreciation of what OP is asking.

From the context provided in your question, "better" likely means size and/or speed of compiled code, perhaps to support a small micro-processor.

The CORDIC ( COordinate Rotation DIgital Computer ) algorithm is very suitable for use on smaller uP, and FPGA implementations that have limited mathematical computation capabilities as it computes the sine and cosine of a value using only basic arithmetic (addition, subtraction and shifts). More about CORDIC, and how to use it to produce sine/cosine of an angle are provided here.

There are also several sites that provide algorithm implementation examples. Simple CORDIC is one that includes detailed explanations on how to generate a table that can then be pre-compiled for use on your target device, as well as code to test the output of the following function (which uses fixed point math):

(See documentation of following, and other functions in link)

#define cordic_1K 0x26DD3B6A
#define half_pi 0x6487ED51
#define MUL 1073741824.000000
#define CORDIC_NTAB 32
int cordic_ctab [] = {0x3243F6A8, 0x1DAC6705, 0x0FADBAFC, 0x07F56EA6, 0x03FEAB76, 0x01FFD55B, 
0x00FFFAAA, 0x007FFF55, 0x003FFFEA, 0x001FFFFD, 0x000FFFFF, 0x0007FFFF, 0x0003FFFF, 
0x0001FFFF, 0x0000FFFF, 0x00007FFF, 0x00003FFF, 0x00001FFF, 0x00000FFF, 0x000007FF, 
0x000003FF, 0x000001FF, 0x000000FF, 0x0000007F, 0x0000003F, 0x0000001F, 0x0000000F, 
0x00000008, 0x00000004, 0x00000002, 0x00000001, 0x00000000 };

void cordic(int theta, int *s, int *c, int n)
{
  int k, d, tx, ty, tz;
  int x=cordic_1K,y=0,z=theta;
  n = (n>CORDIC_NTAB) ? CORDIC_NTAB : n;
  for (k=0; k<n; ++k)
  {
    d = z>>31;
    //get sign. for other architectures, you might want to use the more portable version
    //d = z>=0 ? 0 : -1;
    tx = x - (((y>>k) ^ d) - d);
    ty = y + (((x>>k) ^ d) - d);
    tz = z - ((cordic_ctab[k] ^ d) - d);
    x = tx; y = ty; z = tz;
  }  
 *c = x; *s = y;
}

Edit:
I found the documentation for using the examples at the Simple CORDIC site very easy to follow. However, one small thing I ran into was when compiling the file cordic-test.c the error: use of undeclared identifier 'M_PI' occurred. It appears that when executing the compiled gentable.c file (which generates the cordic-test.c file) the line:

#define M_PI 3.1415926535897932384626

although included in its own declarations, was not included in the printf statements used to produce the file cordic-test.c. Once this was remedied, everything worked as advertised.

As documented, the range of data produced generates 1/4 of a complete sine cycle (-π/2 - π/2 ). The following illustration contains a representation of the actual data produced between the light blue dots. The remainder of the sine signal is fabricated via mirroring and transposing the original data section.

Answer 3

Generating an accurate sine function requires an amount of resource (CPU cycles and memory) that is unwarranted in this application. Your aim to generate a "smooth" sine curve is failing to consider the requirements of the application.

• While when you draw the curve you may observe imperfections, when you apply that curve to an LED PWM drive, the human eye will not perceive those imperfections at all.

• Neither is the human eye likely to perceive the difference in brightness between adjacent values in even a 40-step curve, so interpolation is not necessary.

• It will be more efficient in general if you generate a sine function that generates the appropriate PWM drive values directly without floating point. In fact rather than a sine function a scaled raised cosine would be more appropriate, so that an input of zero results in an output of zero, and an input of half the number of values in the cycle results in the maximum value for your PWM drive.

The following function generates a raised cosine curve for an 8-bit FSD PWM from a 16 value (and 16 bytes) lookup generating a 59-step cycle. So it is both memory and performance efficient compared to your 40 step floating point implementation.

#include <stdint.h>

#define LOOKUP_SIZE 16
#define PWM_THETA_MAX (LOOKUP_SIZE * 4 - 4)

uint8_t RaisedCosine8bit( unsigned n )
{
    static const uint8_t lookup[LOOKUP_SIZE] = { 0, 1, 5, 9,
                                                 14, 21, 28, 36,
                                                 46, 56, 67, 78,
                                                 90, 102, 114, 127} ;
    uint8_t s = 0 ;
    n = n % PWM_THETA_MAX ;

    if( n < LOOKUP_SIZE )
    {
        s = lookup[n] ;
    }
    else if( n < LOOKUP_SIZE * 2 - 1 )
    {
        s = 255 - lookup[LOOKUP_SIZE * 2 - n - 2] ;
    }
    else if( n < LOOKUP_SIZE * 3 - 2 )
    {
        s = 255 - lookup[n - LOOKUP_SIZE * 2 + 2] ;
    }
    else
    {
        s = lookup[LOOKUP_SIZE * 4 - n - 4] ;
    }

    return s ;
}

For an input of 0 <= theta < PWM_THETA_MAX the curve looks like this:

Which is I suggest plenty smooth enough for illumination.

In practice you might use it thus:

for(;;)
{
    for( unsigned i = 0; i < PWM_THETA_MAX; i++ )
    {
        LedPwmDrive( RaisedCosine8bit( i ) ) ;
        Delay( LED_UPDATE_DLEAY ) ;
    }
}

If your PWM range is not 0 to 255, simply scale the output of the function; 8-bit resolution is more than enough for the task.

Answer 4

The classic hack to draw a circle (and hence generate a sine wave too) is Hakmem #149 by Marvin Minsky. E.g.,:

#include <stdio.h>

int main(void)
{
    float x = 1, y = 0;

    const float e = .04;

    for (int i = 0; i < 100; ++i)
    {
        x -= e*y;
        y += e*x;
        printf("%g\n", y);
    }
}

It will be slightly eccentric, not a perfect circle, and you may get some values slightly over 1, but you could adjust by dividing by the maximum or rounding. Also, integer arithmetic can be used, and you can eliminate multiplication/division by using a negative power of two for e, so shift can be used instead.

Answer 5

Have you considered modelling the portion of the sine curve from [0..PI] as a parabola? If the brightness of the LED is only intended to be observed by a human eye, the shapes of the curves ought to be similar enough so that little difference would be detected.

You would just need to figure out the appropriate equation to describe it.

Hmmm, ...

Vertex at (PI/2, 1)

X-axis intersections at (0, 0) and (PI, 0)

f(x) = 1 - K * (x - PI/2) * (x - PI/2)

Where K would be ...

K = 4 / (PI * PI)

Answer 6

For a LED, you could probably just do with 16 or so steps without even interpolating. That said, I can see at least two odd things in your sin1() function:

1) You have 40 data points in sine_table, but you're taking the index x1 modulo 41 of the input. That doesn't seem the right way to handle the periodicity, and lets x1 point past the last index of the array.
2) You're then adding +1, so x2 can be even more over the limits of the array.
3) You're using i in the function, but it's only set in the main program. I can't tell what it's supposed to do, but using a global like that in a simple calculation function seems dirty at minimum. Maybe it's supposed to provide the fractional part for the interpolation, but shouldn't you use phase for that.

Here's a simple interpolator, which seems to work. Adjust to taste.

#include <assert.h>

int A[4] = {100, 200, 400, 800};    
int interpolate(float x)
{
    if (x == 3.00) {
        return A[3];
    }
    if (x > 3) {
        return interpolate(6 - x);
    }
    assert(x >= 0 && x < 3);
    int i = x;
    float frac = x - i;
    return A[i] + frac * (A[i+1] - A[i]);
}

Some arbitrary sample outputs:

interpolate(0.000000) = 100
interpolate(0.250000) = 125
interpolate(0.500000) = 150
interpolate(1.000000) = 200
interpolate(1.500000) = 300
interpolate(2.250000) = 500
interpolate(2.999900) = 799
interpolate(3.000000) = 800
interpolate(3.750000) = 500

_{(I'll leave it to the interested reader to replace all occurrences of 3 with a properly defined symbolic constant, to generalize the function further, and to implement calculating the negative phase too.)}

Answer 7

I would go with Bhaskara I's approximation of a sine function. Using degrees, from 0 to 180, you can approximate the value like so

float Sine0to180(float phase)
{
    return (4.0f * phase) * (180.0f - phase) / (40500.0f - phase * (180.0f - phase));
}

if you want to account for any angle, you'd add

float sine(float phase)
{
    float FactorFor180to360 = -1 * (((int) phase / 180) % 2 );
    float AbsoluteSineValue = Sine0to180(phase - (float)(180 * (int)(phase/180)));
    return AbsoluteSineValue * FactorFor180to360;
}

If you want to do it in radians, you'd add

float SineRads(float phase)
{
    return Sine(phase * 180.0f / 3.1416);
}

Here is a graph showing the points calculated with this approximation and also points calculated with the sine function. You can barely see the approximation points peeking out from under the actual sine points.

Answer 8

Unless your application calls for real precision, don't kill yourself coming up with an algorithm for a 40 point sine or cosine wave. Also, the values in your table should match the range of your LED's pwm input.

That said, I took a look at your code and it's output and figured you weren't interpolating between points. With a little modification, I fixed it, and the error between a excel's sign function and yours is off by a max of about 0.0032 or thereabouts. The change is pretty easy to implement and has been tested using tcc, my personal go-to for C algorithm testing.

First off, I added one more point to your sine array. The last point is set to the same value as the first element in the sine array. This fixes the math in your sine function, in particular when you set x1 to (int)phase%40, and x2 to x1+1. Adding the extra point isn't necessary, as you could set x2 to (x1+1)%40, but I chose the first approach. I'm just pointing out different ways you could accomplish this. I also added the calculation of a remainder (Basically phase - (int)phase). I'm using the remainder for the interpolation. I also added a temporary sine value holder and a delta variable.

const int sine_table[41] = 
{0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767,  32364, 31163, 29196, 26509, 23170, 
19260, 14876, 10125, 5125, 0, -5126, -10126,-14877,
-19261, -23171, -26510, -29197, -31164, -32365, -32768, -32365,
-31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126, 0};

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase)
{
    int tsv,delta;
    float rem;

    rem = phase - (int)phase;
    x1 = (int) phase % 40;
    x2 = (x1 + 1);

    tsv=sine_table[x1];
    delta=sine_table[x2]-tsv;

    y = tsv + (int)(rem*delta);
    return y;
}

int main()
{
    int i;  
    for(i=0;i<420;i++)
    {
       printf("%.2f, %f\n",0.1*i,sin1(0.1*i)/32768);
    }
    return 0;
}

The results look pretty good. Comparing the linear approximation vs the system's floating point sine function gave me the error plot shown below.

Answer 9

You could use the first few terms of the Taylor series expansion of sin. You can use as few terms as needed to reach your intended level of precision - a few more terms than the below example should start to bump up against the limits of a 32-bit float.

Example:

#include <stdio.h>

// Please use the built-in floor function if you can. 
float my_floor(float f) {
    return (float) (int) f;
}

// Please use the built-in fmod function if you can.
float my_fmod(float f, float n) {
    return f - n * my_floor(f / n);
}

// t should be in given in radians.
float sin_t(float t) {
    const float PI = 3.14159265359f;

    // First we clamp t to the interval [0, 2*pi) 
    // because this approximation loses precision for 
    // values of t not close to 0. We do this by 
    // taking fmod(t, 2*pi) because sin is a periodic
    // function with period 2*pi.
    t = my_fmod(t, 2.0f * PI);

    // Next we clamp to [-pi, pi] to get our t as
    // close to 0 as possible. We "reflect" any values 
    // greater than pi by subtracting them from pi. This 
    // works because sin is an odd function and so 
    // sin(-t) = -sin(t), and the particular shape of sin
    // combined with the choice of pi as the endpoint
    // takes care of the negative.
    if (t >= PI) {
        t = PI - t;
    }

    // You can precompute these if you want, but
    // the compiler will probably optimize them out.
    // These are the reciprocals of odd factorials.
    // (1/n! for odd n)
    const float c0 = 1.0f;
    const float c1 = c0 / (2.0f * 3.0f);
    const float c2 = c1 / (4.0f * 5.0f);
    const float c3 = c2 / (6.0f * 7.0f);
    const float c4 = c3 / (8.0f * 9.0f);
    const float c5 = c4 / (10.0f * 11.0f);
    const float c6 = c5 / (12.0f * 13.0f);
    const float c7 = c6 / (14.0f * 15.0f);
    const float c8 = c7 / (16.0f * 17.0f);

    // Increasing odd powers of t.
    const float t3  = t * t * t;
    const float t5  = t3 * t * t;
    const float t7  = t5 * t * t;
    const float t9  = t7 * t * t;
    const float t11 = t9 * t * t;
    const float t13 = t9 * t * t;
    const float t15 = t9 * t * t;
    const float t17 = t9 * t * t;

    return c0 * t - c1 * t3 + c2 * t5 - c3 * t7 + c4 * t9 - c5 * t11 + c6 * t13 - c7 * t15 + c8 * t17;
}

// Test the output
int main() {
    const float PI = 3.14159265359f;
    float t;

    for (t = 0.0f; t < 12.0f * PI; t += (PI * 0.25f)) {
        printf("sin(%f) = %f\n", t, sin_t(t));
    }

    return 0;
}

Example output:

sin(0.000000) = 0.000000
sin(0.785398) = 0.707107
sin(1.570796) = 1.000000
sin(2.356194) = 0.707098
sin(3.141593) = 0.000000
sin(3.926991) = -0.707107
sin(4.712389) = -1.000000
sin(5.497787) = -0.707098
sin(6.283185) = 0.000398
...
sin(31.415936) = 0.000008
sin(32.201332) = 0.707111
sin(32.986729) = 1.000000
sin(33.772125) = 0.707096
sin(34.557522) = -0.000001
sin(35.342918) = -0.707106
sin(36.128315) = -1.000000
sin(36.913712) = -0.707100
sin(37.699108) = 0.000393

As you can see there is still some room for improvement in precision. I am not a genius with floating point arithmetic so probably some of it has to do with the floor/fmod implementations or the specific order the mathematical operations are performed in.

Answer 10

Since you try to generate a signal, i think use a differential equation should not be a bad idea ! it give something like that

#include <stdlib.h>
#include <stdio.h>

#define DT (0.01f) //1/s
#define W0 (3)     //rad/s

int main(void) {
    float a = 0.0f;
    float b = DT * W0;
    float tmp;

    for (int i = 0; i < 400; i++) {
        tmp = (1 / (1 + (DT * DT * W0 * W0))) * (2 * a - b);
        b = a;
        a = tmp;

        printf("%f\n", tmp);
    }
}

Still set the amplitude and frequency of the signal is a pain in the neck :/

Answer 11

It would help if you explain why you don't want the built-in function, but as others have said, the Taylor series is one way to estimate the value. However, the other answers seem to actually be using the Maclaurin series, not Taylor. You should have a lookup table of both sine and cosine. Then find x₀, the closest x value in your lookup table to the x you want, and find d = x-x₀. Then

sin(x) =sin(x₀)+cos(x₀)*d-sin(x₀)*d²/2-cos(x₀)*d³/6 + ...

If your lookup table is such that d<.01, then you'll be getting more than two digits of precision per term.

Another method is to use the fact that if x = x₀+d, then

sin(x) = sin(x₀)*cos(d)+cos(x₀)*sin(d)

You can use a lookup table to get sin(x₀) and cos(x₀), and then use Maclaurin series to get cos(d) and sin(d).

Answer 12

In the case you need float result and accept resolution of whole degrees, sin can be looked up with this code. If you need int (not float) numbers or smaller table, the code can be easily modified:

    /**
     * Sinus lookup table
     * The table only covers whole degrees lookup, and only in I. quarter of the circle. The remaining three quarters
     * of the circle are figured out with logic.
     */
    const float sinus_I_quarter[91] =
    {
        0.0000, 0.0175, 0.0349, 0.0523, 0.0698, 0.0872, 0.1045, 0.1219, 0.1392, 0.1564, // 00 .. 09
        0.1736, 0.1908, 0.2079, 0.2250, 0.2419, 0.2588, 0.2756, 0.2924, 0.3090, 0.3256, // 10 .. 19
        0.3420, 0.3584, 0.3746, 0.3907, 0.4067, 0.4226, 0.4384, 0.4540, 0.4695, 0.4848, // 20 .. 29
        0.5000, 0.5150, 0.5299, 0.5446, 0.5592, 0.5736, 0.5878, 0.6018, 0.6157, 0.6293, // 30 .. 39
        0.6428, 0.6561, 0.6691, 0.6820, 0.6947, 0.7071, 0.7193, 0.7314, 0.7431, 0.7547, // 40 .. 49
        0.7660, 0.7771, 0.7880, 0.7986, 0.8090, 0.8192, 0.8290, 0.8387, 0.8480, 0.8572, // 50 .. 59
        0.8660, 0.8746, 0.8829, 0.8910, 0.8988, 0.9063, 0.9135, 0.9205, 0.9272, 0.9336, // 60 .. 69
        0.9397, 0.9455, 0.9511, 0.9563, 0.9613, 0.9659, 0.9703, 0.9744, 0.9781, 0.9816, // 70 .. 79
        0.9848, 0.9877, 0.9903, 0.9925, 0.9945, 0.9962, 0.9976, 0.9986, 0.9994, 0.9998, // 80 .. 89
        1.0000                                                                          // 90
    };

    #define CIRCLE_QUARTER_1        1
    #define CIRCLE_QUARTER_2        2
    #define CIRCLE_QUARTER_3        3
    #define CIRCLE_QUARTER_4        4
    

    /**
     * The function converts passed angle into 0..89 degrees and looks up a corresponding sin() value.
     * After lookup the return value is adjusted depending what was the original passed angle by figuring
     * out belonging to circle quarter I., II, II. or IV.
     * Any angle outside the accepted angle range of 0 .. 359 will be automatically corrected.
     *
     * @param angle     Whole angle in degrees as used in mathematics: 0 degrees is EAST (in meteorology, 0 degrees wind is NORTH)
     * @return          The looked-up 4-decimal points accurate sinus
     */
    float sinus_lookup (unsigned int angle)
    {
        float sin_value;
        unsigned int circle_quarter;
    
        // correct angles outside the accepted angle range into 0 .. 359
        if (angle > 359u)
            angle = angle % 360u;
    
        circle_quarter = 1 + (angle / 90u);
    
        switch (circle_quarter)
        {
            case CIRCLE_QUARTER_1: // 00 .. 89
                sin_value = sinus_I_quarter[angle];
                break;
    
            case CIRCLE_QUARTER_2: // 90 .. 179
                sin_value = sinus_I_quarter[180 - angle];
                break;
    
            case CIRCLE_QUARTER_3: // 180 .. 269
                sin_value = -sinus_I_quarter[angle - 180];
                break;
    
            case CIRCLE_QUARTER_4: // 270 .. 359
                sin_value = -sinus_I_quarter[360 - angle];
                break;
        }
    
        return sin_value;
    }