Akifumi Takahashi / CubicSpline

This lib is considered to be used as a sensor's calibration program. Calibration with Spline Interpolation might be useful in the case that you want some model expressing relationship such like between a value of physical quantity and your sensor's voltage, but you cannot estimate a model such as liner, square, cubic polynomial, or sine curve. This makes (Parametric) Cubic Spline Polynomial Model (Coefficients of the polynomial) from some sample plots(e.g. sets of (value, voltage)). The inverse function (x,y)->(y,x) has been implemented so as to get analog data (not stepping or leveled data).

Fork of TRP105F_Spline by Akifumi Takahashi

Diff: CubicSpline.cpp

Revision:
12:b3e07a2220bc
Parent:
11:a279e31d8499
Child:
13:9a51747773af
--- a/CubicSpline.cpp	Sun May 29 09:31:44 2016 +0000
+++ b/CubicSpline.cpp	Sun May 29 12:15:19 2016 +0000
@@ -1,4 +1,9 @@
 #define DEBUG
+#define VERSION_C
+#define DEBUG_MAKE_MODEL
+//#define DEBUG_SOLVE
+#define DEBUG_GETX "DEBUG_GETX\n"
+//#define DEBUG_GETY "DEBUG_GETY\n"
 #include "CubicSpline.h"
 
 //  To get voltage of TRP105F
@@ -189,8 +194,6 @@
                    +pow(_Sample_Set[i].y - _Sample_Set[i-1].y, 2));
 }
 
-#define VERSION_C
-#define DEBUG_MAKE_MODEL
 //
 //  Function to define _u_spline, specific constants of spline.
 //
@@ -326,7 +329,6 @@
 //  Function to solve a cubic polinomial
 //  by using Gardano-Tartaglia formula
 //
-#define DEBUG_SOLVE
 void CubicSpline2d::_solve_cubic_f(
     std::complex<double>* arg_t,
     const double  arg_C[4],
@@ -378,8 +380,8 @@
         v = std::complex<double>(-q,-sqrt(-D));
         //u = pow(u, 1/3);
         //v = pow(v, 1/3);
-        u = std::exp(std::log(u)/3);
-        v = std::exp(std::log(u)/3);
+        u = std::exp(std::log(u)/std::complex<double>(3.0,0.0));
+        v = std::exp(std::log(u)/std::complex<double>(3.0,0.0));
     }
     //One real root and two complex root
     else {
@@ -411,7 +413,6 @@
 #endif
 }
 
-#define DEBUG_GETX "DEBUG_GETX\n"
 double  CubicSpline2d::getX(double arg_y)
 {
     double x;
@@ -450,15 +451,15 @@
                 } else*/ if (_Sample_Set[ival].t <= t_sol[i].real() && t_sol[i].real() < _Sample_Set[ival+1].t) {
                     t_real.push_back(t_sol[i].real());
                     t_ival.push_back(ival);
+#ifdef DEBUG_GETX
+                    g_Serial_Signal.printf("(t, i) = (%f, %d)\n", t_real[t_real.size() - 1], ival);
+#endif
                 }
-#ifdef DEBUG_GETX
-                g_Serial_Signal.printf("(t, i) = (%f, %d)\n", t_real[t_real.size() - 1], ival);
-#endif
             }
         }
     }
 
-    if(t_real.size > 0) {
+    if(t_real.size() > 0) {
         the_t = t_real[0];
         the_i = t_ival[0];
         //if t's size is bigger than 1
@@ -469,9 +470,12 @@
             }
         }
     } else {
+#ifdef DEBUG_GETX
+        g_Serial_Signal.printf("LastPoint\n");
+#endif
         the_t = _Last_Point.t;
-        for (int i = 0; i < _Sumple_Num - 1; i++)
-            if(_Sumple_Set[i].t <= the_t && the_t <= _Sample_Set[i+1])
+        for (int i = 0; i < _Sample_Num - 1; i++)
+            if(_Sample_Set[i].t <= the_t && the_t <= _Sample_Set[i+1].t)
                 the_i = i;
     }
     for(int i = 0; i < 4; i++) C[i] = _C_x[i][the_i];
@@ -479,13 +483,15 @@
 #ifdef DEBUG_GETX
     g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)\n",the_t, the_i);
 #endif
-
+    _Last_Point = (Vxyt) {
+        x, arg_y, the_t
+    };
     return x;
 }
 
-#define DEBUG_GETY "DEBUG_GETY\n"
 double  CubicSpline2d::getY(double arg_x)
 {
+
     double y;
     double C[4];
     double the_t;
@@ -506,24 +512,31 @@
     for(int ival = 0; ival < _Sample_Num - 1; ival++) {
         for(int i = 0; i < 4; i++) C[i] = _C_x[i][ival];
         _solve_cubic_f(t_sol, C, arg_x);
+#ifdef DEBUG_GETY
+        g_Serial_Signal.printf("interval:%d \t %f < t < %f\n", ival, _Sample_Set[ival].t, _Sample_Set[ival+1].t);
+#endif
         for(int i = 0; i < 3; i++) {
             //  regarding only real solution
             //  acuracy (error range) is supposed +-10E-3 here(groundless)
             if(std::abs(t_sol[i].imag()) < 0.000001) {
-                /*  */ if (ival == 0               && t_sol[i].real() < _Sample_Set[ival].t) {
+                /*   if (ival == 0               && t_sol[i].real() < _Sample_Set[ival].t) {
                     t_real.push_back(_Sample_Set[ival].t);
                     t_ival.push_back(ival);
                 } else if (ival == _Sample_Num - 2 && _Sample_Set[ival + 1].t <= t_sol[i].real()) {
                     t_real.push_back(_Sample_Set[ival + 1].t);
                     t_ival.push_back(ival);
-                } else if (_Sample_Set[ival].t <= t_sol[i].real() && t_sol[i].real() < _Sample_Set[ival+1].t) {
+                } else*/ if (_Sample_Set[ival].t <= t_sol[i].real() && t_sol[i].real() < _Sample_Set[ival+1].t) {
                     t_real.push_back(t_sol[i].real());
                     t_ival.push_back(ival);
+#ifdef DEBUG_GETY
+                    g_Serial_Signal.printf("(t, i) = (%f, %d)\n", t_real[t_real.size() - 1], ival);
+#endif
                 }
             }
         }
+    }
 
-
+    if(t_real.size() > 0) {
         the_t = t_real[0];
         the_i = t_ival[0];
         //if t's size is bigger than 1
@@ -533,14 +546,23 @@
                 the_i = t_ival[i];
             }
         }
+    } else {
+#ifdef DEBUG_GETY
+        g_Serial_Signal.printf("LastPoint\n");
+#endif
+        the_t = _Last_Point.t;
+        for (int i = 0; i < _Sample_Num - 1; i++)
+            if(_Sample_Set[i].t <= the_t && the_t <= _Sample_Set[i+1].t)
+                the_i = i;
     }
-
     for(int i = 0; i < 4; i++) C[i] = _C_y[i][the_i];
     y = _cubic_f(the_t - _Sample_Set[the_i].t, C);
 #ifdef DEBUG_GETY
-    g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)i\n",the_t, the_i);
+    g_Serial_Signal.printf("(the_t, the_i):= (%f , %d)\n",the_t, the_i);
 #endif
-
+    _Last_Point = (Vxyt) {
+        y, arg_x, the_t
+    };
     return y;
 }
 
@@ -558,6 +580,7 @@
         ft[i]= _Sample_Set[i].x;
     }
     _makeModel(t,ft,_C_x);
+    
     for(int i = 0; i < _Sample_Num; i++) {
         ft[i]= _Sample_Set[i].y;
     }
@@ -739,15 +762,19 @@
 void CubicSpline2d::printOutData()
 {
     FILE *fp;
+    double x,y;
     double d = (_Sample_Set[_Sample_Num - 1].x - _Sample_Set[0].x) / 100;
 
     fp = fopen("/local/log.txt", "w");  // open file in writing mode
 
-    fprintf(fp, "x, y\n");
-    for(int ival = 0; ival < _Sample_Num; ival++) {
-        for(double x = _Sample_Set[ival].x; x < _Sample_Set[ival+1].x; x += d) {
-            fprintf(fp, "%f,%f\n", x, getY(x));
+    fprintf(fp, "x, y, (t)\n");
+    for(int ival = 0; ival < _Sample_Num - 1; ival++) {
+        fprintf(fp, "ival: %d \n", ival);
+        for(x = _Sample_Set[ival].x; x < _Sample_Set[ival + 1].x; x += d) {
+            y = getY(x);
+            fprintf(fp, "%f,%f,(%f)\n", x, y, sqrt(x*x + y*y));
         }
+        fprintf(fp, "ival: %d \n", ival);
     }
 
     fprintf(fp, "\nSample:dst, vol\n");