Mario.Tapilouw

Saturday, March 03, 2012

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!

Labels: , , ,

9 Comments:

  • Dear Mario,

    I was taking a look to the algorithm you posted since I need to fit a sinusoidal curve to a set of data.

    Then I noticed that the code corresponding to the entry labeled "Then fill the data with a sine or cosine function and the simulated input data:" is incomplete.

    Would you be so kind as to fill the missing parts?

    And thank you very much for the algorithm.

    Regards,

    José.

    By Blogger Jose Navarro, at April 25, 2018 at 5:53 PM  

  • Dear Mario,

    first of all, thanks for the algorithm you have made publicly available. I hope to use it in a problem related to sinusoidal curve fitting I have.

    However, there is a problem: the code related to the part where you wrote "Then fill the data with a sine or cosine function and the simulated input data:" is missing.

    Could you please fill this part so the algorithm is complete?

    Thank you very much!

    José.

    By Blogger Jose Navarro, at April 25, 2018 at 5:55 PM  

  • Hi José,

    Thanks for stopping by. Sorry I didn't realize that it was missing, I'll search my repositories first.

    By Blogger mario, at April 26, 2018 at 9:02 AM  

  • Hi

    Here it is, hope it's not deleted again..

    for(int i=0;i<numOfData;i++)
    {
    x[i] = cos((2*M_PI) * float(i%numOfData) / float(numOfData));
    y[i] = sin((2*M_PI) * float(i%numOfData) / float(numOfData));
    z[i] = your_data[i];
    }

    By Blogger mario, at April 26, 2018 at 9:54 AM  

  • Hi!

    Thank you very much for posting the missing code. Just one final remark: it is shown in your answer to my comment, but it doesn't show up in the webpage itself (I mean, in the area where you explain, step by step, the different actions to take).

    Thanks again,

    José

    By Blogger Jose Navarro, at April 26, 2018 at 4:15 PM  

  • Hi, Mario, it's José once again.

    I'm afraid that there's another snippet of code missing. This time, in the box where you write the code labeled "Prepare the variables for performing least square fitting.". At least, it looks like it's missing. I have included BELOW the code I have copied from the webpage so you can check the problem by yourself.

    One last question: once this issue is clarified, may I say that just by contatenating the several snippets of code I'll have the WHOLE algorithm to fit a sinusoidal to my data series?

    Thanks,

    José.

    -----------

    // 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 // <------------------------ MISSING CODE AGAIN

    By Blogger Jose Navarro, at April 26, 2018 at 4:39 PM  

  • Fixed that missing part, thanks José

    By Blogger mario, at April 26, 2018 at 5:03 PM  

  • Mario,

    thanks for posting the missing parts of the code!

    Just one last question: I have checked the code; to do this, I have created a synthetic sinusoidal wave using a function of the type y = A + C * sin(b*x), that is, the same kind of function you start describing in the post.

    Then I run the program and get the magnitude (I assume that this is the AMPLITUDE, that is, C in the formula above) and the phase, that is, b in the same formula.

    The values of C and b do not look like the original ones used to generate the synthetic data (C goes up to the sky, being something like 1.5e+14 while the original value was just 5.). And I don't know where to obtain A, the constant term.

    The only thing I did was to generate the synthetic sinusoidal wave populating vector dataz and then stating the number of elements in this array by means of variable numOfData.

    Is this the right way to proceed? What am I doing wrong?

    Thank you very much.

    José.

    By Blogger Jose Navarro, at April 26, 2018 at 10:29 PM  

  • I don't know what's wrong with your code, do you mind sending the source code to me? mario.tapilouw (at) gmail (dot) com

    By Blogger mario, at April 27, 2018 at 8:57 AM  

Post a Comment

<< Home