Sine Fitting

Sine is not a linear function, but fitting series of data to sine function is actually not a difficult task. I need to apply this fitting function so I was trying to find a way to do this.

Sine function can be represented in general form as:

- y(x) = A + C * sin(x + b)

this function can be rewritten as:

- y(x) = A + C * sin(b) * cos(x) + C * cos(b) * sin(x)

or can also be written as:

- y = b0 + b1 * x1 + b2 * x2

which is a general linear regression form and can be solved by using LM function to obtain b0, b1, and b2, then based on these values we can calculate A, b and C

1. A = b0

2. b1^2 + b2^2 = c^2 * (sin^2(b) + cos^2(b)) = c^2

C = sqrt(b1^2 + b2^2)

3. tan(b) = b1/b2

b = arctan(b1/b2)

C++ implementation:

1. Prepare the buffer for datax, datay for for cos(b), sin(b), and dataz for the input data ) and variables for the result, phase and mag:

float * datax = new float[numOfData]; 
float * datay = new float[numOfData]; 
float * dataz = new float[numOfData];  
float phase(0); 
float mag(0);

2. Then fill the data with a sine or cosine function and the simulated input data:

for(int i=0;i < numOfData;i++)
{
  datax[i] = cos((2*M_PI) * float(i%numOfData) / float(numOfData));
  datay[i] = sin((2*M_PI) * float(i%numOfData) / float(numOfData));
  dataz[i] = yourData[i];
}

3. Prepare the variables for performing least square fitting.

float a(0); 
float b(0); 
float c(0); 
float p[3] = {0}; // product of fitting equation 
float XiYi(0); 
float XiZi(0); 
float YiZi(0); 
float XiXi(0); 
float YiYi(0); 
float Xi(0); 
float Yi(0); 
float Zi(0);  

for(int i=0; i < numOfData; i++)
  {
    XiYi += datax[i] * datay[i];
    XiZi += datax[i] * dataz[i];
    YiZi += datay[i] * dataz[i];
    XiXi += datax[i] * datax[i];
    YiYi += datax[i] * datay[i];
    Xi   += datax[i];
    Yi   += datay[i];
    Zi   += dataz[i];
  }


4. Then calculate the inverse of the matrix as follows:

float A[3][3];
float B[3][3];
float C[3][3];
float X[3][3];
int i;
int j;
float x = 0;
float n = 0; //n is the determinant of A

A[0][0] = XiXi;
A[0][1] = XiYi;
A[0][2] = Xi;
A[1][0] = XiYi;
A[1][1] = YiYi;
A[1][2] = Yi;
A[2][0] = Xi;
A[2][1] = Yi;
A[2][2] = numOfData;

for( i=0;i<3;i++)
{
 for(int j=0;j<3;j++)
 {
   B[i][j] = 0;
   C[i][j] = 0;
 }
}

for( i=0, j=0;j<3;j++)
{
 if(j == 2)
 {
   n += A[i][j] * A[i+1][0] * A[i+2][1];
 }
 else if(j == 1)
 {
   n += A[i][j] * A[i+1][j+1] * A[i+2][0];
 }
 else
 {
   n += A[i][j] * A[i+1][j+1] * A[i+2][j+2];
 }
}

for( i=2, j=0;j<3;j++)
{
 if(j == 2)
 {
   n -= A[i][j] * A[i-1][0] * A[i-2][1];
 }
 else if(j == 1)
 {
   n -= A[i][j] * A[i-1][j+1] * A[i-2][0];
 }
 else
 {
   n -= A[i][j]*A[i-1][j+1]*A[i-2][j+2];
 }
}

// Check determinant n of matrix A
if (n)
{
 x = 1.0/n;

 for(i=0;i<3;i++)
 {
   for(j=0;j<3;j++)
   {
     B[i][j] = A[j][i];
   }
 }

 C[0][0] =  (B[1][1] * B[2][2]) - (B[2][1] * B[1][2]);
 C[0][1] = ((B[1][0] * B[2][2]) - (B[2][0] * B[1][2])) * (-1);
 C[0][2] =  (B[1][0] * B[2][1]) - (B[2][0] * B[1][1]);

 C[1][0] = ((B[0][1] * B[2][2]) - (B[2][1] * B[0][2])) * (-1);
 C[1][1] =  (B[0][0] * B[2][2]) - (B[2][0] * B[0][2]);
 C[1][2] = ((B[0][0] * B[2][1]) - (B[2][0] * B[0][1])) * (-1);

 C[2][0] =  (B[0][1] * B[1][2]) - (B[1][1] * B[0][2]);
 C[2][1] = ((B[0][0] * B[1][2]) - (B[1][0] * B[0][2])) * (-1);
 C[2][2] =  (B[0][0] * B[1][1]) - (B[1][0] * B[0][1]);

 for( i=0;i<3;i++)
 {
   for( j=0;j<3;j++)
   {
     X[i][j] = C[i][j] * x;
   }
 }

 p[0] = XiZi;
 p[1] = YiZi;
 p[2] = Zi;

 a = X[0][0] * p[0] + X[0][1] * p[1] + X[0][2] * p[2];
 b = X[1][0] * p[0] + X[1][1] * p[1] + X[1][2] * p[2];
 c = X[2][0] * p[0] + X[2][1] * p[1] + X[2][2] * p[2];
}
else  // determinant=0
{
 a = 1;
 b = 1;
 c = 0; 
}

5. Then we could obtain the magnitude:

mag = sqrt(a*a + b*b);

6. Finally we could calculate the phase of the signal:

float phi;

if((b == 0) && (a >= 0))
{
 phase = M_PI/2;
}
else if((b == 0) && (a >= 0))
{
 phase = 3*M_PI/2;
}
else if ((a >= 0) && (b > 0))
{
 phi = atan(fabs(a/b));
 phase = phi;
}
else if ((a >= 0) && (b < 0))
{
 phi = atan(fabs(a/b));
 phase = M_PI - phi;
}
else if ((a <  0) && (b > 0))
{
 phi = atan(fabs(a/b));
 phase = 2 * M_PI - phi;
}
else if ((a <  0) && (b < 0))
{
 phi = atan(fabs(a/b)); 
 phase = phi + M_PI;
}

Then it's time to visualize the result...

Good luck!

GNU Scientific Library

Recently, I spend some time in looking for libraries that can fasten my software development. The reason is simple, if I want to develop my software for example doing a common algorithm, I don't need to spend time for developing the algorithm which might not be optimized and unstable.

The fact is that many people out there who share the same need and they have already created a group for developing an open source softwares/libraries. The advantage of such a group is that many people evaluates the software and many people will criticize when there's a bug in the software while the others will contribute to resolve the bug. Quite wonderful right? You have a lot of bright people from around the world assisting you (not exaggerating but it's quite true..).

GNU Scientific Library provides some scientific calculation such as linear least squares, FFT, Linear Algebra, Numerical Differentiation, etc. You might want to see this website for the available routines provided by GSL. It helps you focus on your task and appreciate others who has contributed their time developing this library.

After spending the whole afternoon configuring my compiler to work with this library, got a lot of error messages which are discouraging, finally I found the library that has been ported into Windows platform here

Before using this library, the compiler has to be configured. The configuration includes setting the additional include path and the library path. After that your compiler is ready to use. My compiler is Borland C++ Builder, you need to convert the library into the Borland format instead of using Ms VC++ format (see my note about implib)

Finally, I tried this sample code from the documentation:

#include < stdio.h >

#include < gsl/gsl_sf_bessel.h >

int main (void)
{
double x = 5.0;
double y = gsl_sf_bessel_J0 (x);
printf ("J0(%g) = %.18e\n", x, y);
return 0;
}

and it works :D

I plan to use this library for my projects, it will reduce my work in implementing the basic functions.