Optimization of harmonics calculation

Question

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.

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. — Cedron Dawg, Jul 21 '20 at 13:33

Cedron Dawg · Accepted Answer · 2020-07-20T23:59:14.427

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:

Angle Addition Formulas from Euler's Formula

Optimization of harmonics calculation

1 Answers1

Linked