# Optimization of harmonics calculation

I need to compute the sine and cosine of an argument along with n "harmonics"

$$\begin{matrix} \sin(x) & \cos(x) \\ \sin(2x) & \cos(2x) \\ \cdots \\ \sin(nx) & \cos(nx) \end{matrix}$$

This takes a lot of computing time and since this is a real-time system, I need to optimize this calculation.

Here's my solution so far :

I defined a rotation matrix $$R = \begin{pmatrix} \cos(x) \ -\sin(x) \\ \sin(x) \ \cos(x) \\ \end{pmatrix}$$

I compute $$[\cos(nx); \sin(nx)] = R [\cos((n-1)x); \sin((n-1)x)]$$; I simply iterate to compute the missing $$n-1$$ harmonics.

This method is 3 times fast than computing all harmonics separately.

I'm curious, are there better solutions? I know a look-up table could be even faster but I would prefer not to go this way.

• I'll try to improve the math formating later.
– Ben
Jul 20 '20 at 14:45
• For the kids playing along at home. You can create your very own trig table by starting with one degree ( $2\pi/360$ radians ). Use Taylor series to get the initial values, then apply this method up to 45 degrees. Replicate from there. ;-) Extra credit: Use double angle formulas and known values to minimize error accumulation. Jul 21 '20 at 13:33

No, you can't do it much simpler or faster than that. In the old days, subscripts used to cost time, not sure that is true any more.

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

//=========================================================
int main( int argCount, char *argValues[] )
{
double x = 1.2;

double C = cos( x );
double S = sin( x );

double CA[10];
double SA[10];

CA[0] = 1.0;
SA[0] = 0.0;

CA[1] = C;
SA[1] = S;

for( int h = 2; h < 10; h++ )
{
CA[h] = CA[h-1] * C - SA[h-1] * S;
SA[h] = CA[h-1] * S + SA[h-1] * C;

printf( "%2d   %10.6f %10.6f      %10.6f %10.6f\n",
h, CA[h], SA[h], cos( h * x ), sin( h * x ) );
}

return 0;
}
//=========================================================


Results:

 2    -0.737394   0.675463       -0.737394   0.675463
3    -0.896758  -0.442520       -0.896758  -0.442520
4     0.087499  -0.996165        0.087499  -0.996165
5     0.960170  -0.279415        0.960170  -0.279415
6     0.608351   0.793668        0.608351   0.793668
7    -0.519289   0.854599       -0.519289   0.854599
8    -0.984688  -0.174327       -0.984688  -0.174327
9    -0.194330  -0.980936       -0.194330  -0.980936


This implementation is based on complex multiplication rather than matrix rotation. In this case, conceptually they are different, mechanically they are identical.

You will accumulate a little bit of error, but for small sets this remains in the $$10^{-17}$$ region.

See: