### 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`

`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`

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: least square, math, noisy data, sine fitting