Function for approximating a sine signal

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.

• Any specific microcontroller? Dec 22, 2017 at 20:19
• The phase is given as thousandths of a radian? Dec 22, 2017 at 20:21
• Looks like 1000 phase = 1 revolution, so a phase is 0.36 degrees, or 18/pi radians. Dec 22, 2017 at 20:32
• is phase guaranteed to be positive? Dec 22, 2017 at 22:27
• Just as a general rule, please indent... Dec 23, 2017 at 17:51

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 [0RR] 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;
}

• Thank you very much! That will surely help me a lot. I'm just trying to understand and compile your code. Could you please explain what exactly "unsigned_mult16x16_into_32" means and how to implement this? Dec 23, 2017 at 12:22
• @Peter123 Judging by the name, it multiplies two uint_fast16_ts and return a uint_fast32_t. You can implement it as such, anyway. Dec 23, 2017 at 12:52
• Peter123 @wizzwizz4 is correct. unsigned_mult16x16_into_32() is a fictional name for a function that multiples 2 unsigned 16-bit integers and forms the 32 bit unsigned product. Embedded compilers often have such an optimal function. Else could use (uint32_t)a*b or the like. Dec 23, 2017 at 17:00
• Thanks! In my application, I sometimes need to get an output range mapped to [-a a] with an arbitrary positive integer a. Do you see any possibility to do that without using floating point math? Dec 23, 2017 at 17:17
• @Peter123 Let P = 1024, the period and int_fast16_t sn = ... / m;, with m==4, the amplitude a is 16384. For a different amplitude, m = (P/4)*(P/4)/a; If quotient m is not a whole number, use some integer scaling. Dec 23, 2017 at 19:19

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.

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?

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).

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)

• BTW compared to a lookup table, you would still need to multiply by amplitude, so I don't bother with a LUT. I just calculate the sine.
– RAB
Dec 24, 2017 at 5:19
• PS this sine wave may start to diverge after many cycles due to numerical error so it's best to run it for one cycle only then re-initialize.
– RAB
Dec 24, 2017 at 5:21