Function for approximating a sine signal

Question

I want to approximate a sine signal with a parabola, which is given by the function:

$$f(x) = 0.4053 \cdot x \cdot (3.1415 - x) \space\text{for}\space 0 \le x \le 3.1415$$

In a first approach, I wrote the following code which works fine:

float sin1(int phase)
{
    if(phase % 1000 <= 500)
    {
    return(0.4053*2*3.1415*(phase % 1000)/1000*(3.1415 - 2*3.1415*(phase % 1000)/1000));
    }
    else
    {
    return(-0.4053*2*3.1415*(phase % 500)/1000*(3.1415 - 2*3.1415*(phase % 500)/1000));
    }
}

As I want to use this function on a microcontroller, I'm interested in how to optimize this code regarding memory consumption and performance.

Answer 1

how to optimize this code regarding memory consumption and performance. (?)

1. Consider redefining the angle measurement. It appears int phase is 0.36 degrees. A natural choice would be 1024 phase to 1 revolution. Now the problem becomes one using a binary angle measure. Not only does this simplify small things like scaling, it opens up coding choices as there are many trig function code bases designed withs BAMS, Binary Angular Measurement System.

2. Code uses double constants and so makes for likely slower computation on an embedded system. If all that is needed is float, uses float constants and float math.

   // return(0.4053*2*3.1415*(phase % 1000)/1000*(3.1415 - ...
   return(0.4053f*2*3.1415f*(phase % 1000)/1000*(3.1415f ...
   

3. OP reports phase guaranteed to be positive. yet it is minor to adjust code for full int range.

   // if(phase % 1000 <= 500)
   phase %= 1000;
   if (phase < 0) phase += 1000;
   if (phase <= 500)
   

4. Accuracy: Using 3.1415 for machine_pi can be improved at no computational cost by using #define MACHINE_PI 3.141592653589793f. This will improve the end points in the 5th or 6th decimal place. Use 0.405284734572655f instead of 0.4053 to insure sine(90° (250 phase)) is 1.0f.

5. If your processor lacks a built-in FP hardware, I recommend a fast integer only solution. e.g. CORDIC or a small linear interpolation should only low precision be needed.

---

Using a parabola model like OP's (1 cycle = 1024 BAM). The below code is effectively the same preciseness as OP's code with only integer math: and, or, +, -, *, negate and shift.

Input range [0^R 2π^R] mapped to [0 to 1024]
Output range [-1.0 1.0] mapped to [-16,384 16,384]

int_fast16_t sin_bam(int angle_bam) {
  angle_bam %= 1024;
  if (angle_bam < 0)
    angle_bam += 1024;
  uint_fast16_t sub_angle = angle_bam % (1024u / 2);
  int_fast16_t sn = unsigned_mult16x16_into_32(sub_angle, (1024/2 - sub_angle)) / 4;
  if (angle_bam & 512) {
    sn = -sn;
  }
  return sn;
}

Answer 2

here are the programming principles you ought to be applying:

1. Don't repeat yourself

   Note that the two statements (statement-true & statement-false) are just about the same as each other. If there's a mistake in one and you fix it, you might not fix it correctly in the other. So you really should have just one copy of the formula. This is even more critical on a microcontroller as you want to have the shortest code, right?

2. Don't repeat yourself

   Hmm, sounds familiar. Looks like you have phase % 1000 multiple places. You should do it just once, into a temporary variable e.g. int phase_mod_1000, unless you are completely critical about avoiding temporary storage. I bet that freeing up that float will be worth it. And also note, phase % 500 is just going to be phase_mod_1000 - 500. So three divisions turn into one. The optimizer might be able to do this automatically; this is really just to simplify the code. You should probably have constants (#define) for PHASES_PER_CYCLE and PHASES_PER_HALFCYCLE.

3. Don't do additional work

   I'm sure someone with a good sense of algebra could simplify the arithmetic, but guess what, the optimizer can too. Your compiler has a good optimizer, right?

Answer 3

I'm interested in how to optimize this code regarding memory consumption and performance.

Well, which is it? Are you optimizing for speed, ram, or executable size? You say this is for an embedded project, so it’s likely you’ll need to pick one.

Once I needed to calculate a sine wave for a motion controller. The timing was critical, but we had plenty of ROM and enough RAM that we just included a pre-calculated table of sin values. They never change, so there’s no need to perform an expensive calculation at runtime, particularly if you’re in a timing sensitive function.

I honestly can’t tell if you’re calculating the same values or not, but the principle still applies I think. You only need this code long enough to calculate the values once, so it’s quality is a moot point (so long as it’s bug free).

Answer 4

If you want a sine there is a better way with only one multiply per step.

Let const P = the phase increment = 2 * pi * f * dt
Let const A = the amplitude
Let const H = 2 cos(P) (you only need to calculate H once, maybe by hand)

(note H is slightly less than 2)

Initialize V=0 and W=A*cos(P) = A*H/2

Loop: output V, X = W*H-V, V=W, W=X

So the only multiply is W*H, plus one subtraction and three assignments.

V and W are your first and second points, X is the next point etc. You do have to take a cosine once, to initialize H

If you need a lot of cosines you might want to check out the CORDIC algorithm (see Wikipedia)