fit a 6-th degree curve to sin for coefficients
This commit is contained in:
parent
873e8c4cb7
commit
3cab32b171
30
fit_curve.py
Normal file
30
fit_curve.py
Normal file
@ -0,0 +1,30 @@
|
|||||||
|
import numpy as np
|
||||||
|
from scipy.optimize import curve_fit
|
||||||
|
from matplotlib import pyplot as plt
|
||||||
|
import math
|
||||||
|
|
||||||
|
x = np.linspace(0, math.pi / 2, num=1000)
|
||||||
|
y = np.sin(x)
|
||||||
|
|
||||||
|
|
||||||
|
def test(x, a, b, c, d, e, f):
|
||||||
|
return a * x**6 + b * x**5 + c * x**4 + d * x**3 + e * x**2 + f * x
|
||||||
|
|
||||||
|
|
||||||
|
param, param_cov = curve_fit(test, x, y)
|
||||||
|
|
||||||
|
print("Sine function coefficients:")
|
||||||
|
print(*param)
|
||||||
|
|
||||||
|
|
||||||
|
# ans stores the new y-data according to
|
||||||
|
# the coefficients given by curve-fit() function
|
||||||
|
ans = test(x, *param)
|
||||||
|
|
||||||
|
print("Max error:")
|
||||||
|
print(abs(y-ans).max())
|
||||||
|
|
||||||
|
plt.plot(x, y, 'o', color='red', label="data")
|
||||||
|
plt.plot(x, ans, '--', color='blue', label="optimized data")
|
||||||
|
plt.legend()
|
||||||
|
plt.show()
|
@ -16,14 +16,78 @@ struct math_spec {
|
|||||||
f64 to;
|
f64 to;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
static f64 square(f64 num)
|
||||||
|
{
|
||||||
|
return num * num;
|
||||||
|
}
|
||||||
|
|
||||||
static f64 my_sin(f64 num)
|
static f64 my_sin(f64 num)
|
||||||
{
|
{
|
||||||
return num;
|
f64 sign = 1;
|
||||||
|
if (num < 0) {
|
||||||
|
sign = -1;
|
||||||
|
num = -num;
|
||||||
|
}
|
||||||
|
|
||||||
|
f64 result = (4 / M_PI) * num - square((2 / M_PI) * num);
|
||||||
|
return result * sign;
|
||||||
|
}
|
||||||
|
|
||||||
|
static f64 my_sin_quarter(f64 num)
|
||||||
|
{
|
||||||
|
f64 sign = 1;
|
||||||
|
if (num < 0) {
|
||||||
|
sign *= -1;
|
||||||
|
num = -num;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (num > M_PI / 2) {
|
||||||
|
num = M_PI - num;
|
||||||
|
}
|
||||||
|
|
||||||
|
f64 result = (4 / M_PI) * num - square((2 / M_PI) * num);
|
||||||
|
return result * sign;
|
||||||
|
}
|
||||||
|
|
||||||
|
static f64 my_sin_quarter_fitted(f64 num)
|
||||||
|
{
|
||||||
|
f64 sign = 1;
|
||||||
|
if (num < 0) {
|
||||||
|
sign *= -1;
|
||||||
|
num = -num;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (num > M_PI / 2) {
|
||||||
|
num = M_PI - num;
|
||||||
|
}
|
||||||
|
|
||||||
|
f64 coeffs[] = {
|
||||||
|
-0.0009462013574375033,
|
||||||
|
0.010206999409602957,
|
||||||
|
-0.0018834448707133052,
|
||||||
|
-0.16568577654109595,
|
||||||
|
-0.00024161815691126504,
|
||||||
|
1.0000207681130697
|
||||||
|
};
|
||||||
|
|
||||||
|
f64 result = 0;
|
||||||
|
result += coeffs[0] * num * num * num * num * num * num;
|
||||||
|
result += coeffs[1] * num * num * num * num * num;
|
||||||
|
result += coeffs[2] * num * num * num * num;
|
||||||
|
result += coeffs[3] * num * num * num;
|
||||||
|
result += coeffs[4] * num * num;
|
||||||
|
result += coeffs[5] * num;
|
||||||
|
return result * sign;
|
||||||
}
|
}
|
||||||
|
|
||||||
static f64 my_cos(f64 num)
|
static f64 my_cos(f64 num)
|
||||||
{
|
{
|
||||||
return num;
|
if (num < 0) {
|
||||||
|
num *= -1;
|
||||||
|
}
|
||||||
|
|
||||||
|
f64 result = 1 - square((2 / M_PI) * num);
|
||||||
|
return result;
|
||||||
}
|
}
|
||||||
|
|
||||||
static f64 my_sqrt1(f64 num)
|
static f64 my_sqrt1(f64 num)
|
||||||
@ -60,6 +124,20 @@ int main()
|
|||||||
.from = -M_PI,
|
.from = -M_PI,
|
||||||
.to = M_PI
|
.to = M_PI
|
||||||
},
|
},
|
||||||
|
{
|
||||||
|
.name = "sin_quater",
|
||||||
|
.mine = my_sin_quarter,
|
||||||
|
.reference = sin,
|
||||||
|
.from = -M_PI,
|
||||||
|
.to = M_PI
|
||||||
|
},
|
||||||
|
{
|
||||||
|
.name = "sin_quater_fit",
|
||||||
|
.mine = my_sin_quarter_fitted,
|
||||||
|
.reference = sin,
|
||||||
|
.from = -M_PI,
|
||||||
|
.to = M_PI
|
||||||
|
},
|
||||||
{
|
{
|
||||||
.name = "cos",
|
.name = "cos",
|
||||||
.mine = my_cos,
|
.mine = my_cos,
|
||||||
@ -94,6 +172,7 @@ int main()
|
|||||||
struct math_spec spec = specs[i];
|
struct math_spec spec = specs[i];
|
||||||
|
|
||||||
f64 max_error = 0;
|
f64 max_error = 0;
|
||||||
|
f64 max_error_at = 0;
|
||||||
|
|
||||||
for (int j = 0; j < samples; j++) {
|
for (int j = 0; j < samples; j++) {
|
||||||
f64 input = ((f64)j / samples) * (spec.to - spec.from) + spec.from;
|
f64 input = ((f64)j / samples) * (spec.to - spec.from) + spec.from;
|
||||||
@ -107,10 +186,11 @@ int main()
|
|||||||
|
|
||||||
if (error > max_error) {
|
if (error > max_error) {
|
||||||
max_error = error;
|
max_error = error;
|
||||||
|
max_error_at = input;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
printf("%s max error %f\n", spec.name, max_error);
|
printf("%20s max error %f at %f\n", spec.name, max_error, max_error_at);
|
||||||
}
|
}
|
||||||
|
|
||||||
return 0;
|
return 0;
|
||||||
|
Loading…
Reference in New Issue
Block a user