Recursively create a sine wave given a single sine wave value and the period

Question

I am trying to write a .oct function for Octave that, given a single sine wave value, between -1 and 1, and sine wave period, returns a sine wave vector of period length with the last value in the vector being the given sine wave value. My code so far is:

#include <octave/oct.h>
#include <octave/dColVector.h>
#include <math.h>
#define PI 3.14159265

DEFUN_DLD (sinewave_recreate, args, , "args(0) sinewave value, args(1) is period")
{
octave_value_list retval;

double sinewave_value = args(0).double_value (); 
double period = args(1).double_value ();  
ColumnVector output_sinewave(period);                
double degrees_inc = 360 / period;
double output_sinewave_degrees;

output_sinewave_degrees = asin( sinewave_value ) * 180 / PI;
output_sinewave(period-1) = sin( output_sinewave_degrees * PI / 180 );

for (octave_idx_type ii (1); ii < period; ii++) // Start the loop
   {
   output_sinewave_degrees = output_sinewave_degrees - degrees_inc;

   if ( output_sinewave_degrees < 0 )
   {
   output_sinewave_degrees += 360 ;
   }  

   output_sinewave( period-1-ii ) = sin( output_sinewave_degrees * PI / 180 );
   }

retval(0) = output_sinewave;                                                          

return retval;                                                                        
}

but is giving patchy results. By this I mean that it sometimes recreates the sine wave quite accurately and other times it is way off. I have determined this simply by creating a given sine wave, taking the last value in time and plugging this into the function to recreate the sine wave backwards through time and then comparing plots of the two. Obviously I am doing something wrong, but I can't seem to identify what.

Answer 1

You might try an easier way to go through. Just recall that if

y = sin(x)

then first derivative of y will be equal to

dy/dx = cos(x)

So at every step of computation you add to the current value of y some delta equal to

dy = cos(x) * dx

But that might cut your accuracy down as a side-effect. You could probe it whatever. HTH.

---

It seems that slightly improved equation tend to be more accurate:

dy = cos(x + dx/2) * dx

Take a look at this.

Answer 2

Lets start with some trigonometric identities:

sin(x)^2 + cos(x)^2 == 1
sin(x+y) == sin(x)*cos(y) + sin(y)*cos(x)
cos(x+y) == cos(x)*cos(y) - sin(x)*sin(y)

Given the sine and cosine at a point x, we can exactly calculate the values after a step of size d, after precalculating sd = sin(d) and cd = cos(d):

sin(x+d) = sin(x)*cd + cos(x)*sd
cos(x+d) = cos(x)*cd - sin(x)*sd

Given the initial sine value, you can calculate the initial cosine value:

cos(x) = sqrt(1 - sin(x)^2)

Note that there are two possible solutions, corresponding to the two possible square-root values. Also note that all the angles in these identities are in radians, and d needs to be negative if you're going back through the wave.

Answer 3

Mike's note that there are two possible solutions for cos(x) made me realise that I would need to resolve the phase ambiguity of the sine wave. My second, successful attempt at this function is:

#include <octave/oct.h>
#include <octave/dColVector.h>
#include <math.h>
#define PI 3.14159265

DEFUN_DLD (sinewave_recreate_3, args, , "args(0) sinewave value, args(1) is period, args(2) is the phase")
{
octave_value_list retval;

double sinewave_value = args(0).double_value (); 
double period = args(1).double_value ();  
double phase = args(2).double_value ();
ColumnVector output_sinewave(period);                
double X0 = asin(sinewave_value);

if (sinewave_value < 0 & phase > 180 & phase < 270)
 {
 X0 = PI + (0 - X0); 
 }

if (sinewave_value < 0 & phase >= 270)
 {
 X0 = X0 + 2 * PI; 
 }

if (sinewave_value > 0 & phase > 90)
 {
 X0 = PI - X0; 
 }

if (sinewave_value > 0 & phase < 0)
 {
 X0 = X0 + PI / 2; 
 }

double dx = PI / 180 * (360/period);

for (octave_idx_type ii (0); ii < period; ii++) // Start the loop
 {
 output_sinewave(period-1-ii) = sin(X0 - dx * ii); 
 }

retval(0) = output_sinewave;                                                          

return retval;                                                                        
}

Thanks are also due to Keynslug.