1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ I'll try to improve the math formating later. $\endgroup$
    – Ben
    Jul 20 '20 at 14:45
  • $\begingroup$ 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. $\endgroup$ Jul 21 '20 at 13:33
1
$\begingroup$

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:

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.