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.

Cedron Dawg · Answer

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.6fn",
                    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

One Answer

Add your own answers!

Ask a Question