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
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: , , ,

0 Comments:

Post a Comment

<< Home