interp.h

来自「小波提升格式的源代码」· C头文件 代码 · 共 199 行

/**
   Support for four point polynomial interpolation, using the Lagrange
   formula.

 */
class interpolation {

   private:
   double coef_minus2[4];
   double coef_minus1[4];
   /** lagrange coefficients for x = -0.5 */
   double coef_minus05[4];
   /** lagrange coefficients for x =  0.5 */
   double coef_05[4];
   /** lagrange coefficients for x =  1.5 */
   double coef_15[4];
   /** lagrange coefficients for x =  2.5 */
   double coef_25[4];
   /** lagrange coefficients for x =  3.5 */
   double coef_35[4];
   double coef_40[4];
   double coef_50[4];

   /**
      Evaluate the Lagrange formula at <tt>x</tt> producing a set of 
      coefficients for the equation
      <pre>
      f(x) = c<sub>0</sub>y<sub>0</sub> + c<sub>1</sub>y<sub>1</sub> + c<sub>2</sub>y<sub>2</sub> + c<sub>2</sub>y<sub>2</sub>;
      </pre>

      Where <tt>f(x)</tt> is the interpolated point and
      {y<sub>0</sub>, y<sub>1</sub>, y<sub>2</sub>, y<sub>3</sub>} are
      four points from a time series.

      The reference I used for this is Chapter 3.1 of Numerical
      Recipes in C.  In Numerical Recipies they criticize the Lagrange
      polynomial method because it does not calculate an error value
      for the calculated point.  For some sets of four points the
      value calculated will have a high error and is unreliable.  But
      for the purposes of wavelets the error is not really useful.

      Calculating the coefficients in advance for a four point
      Lagrange polynomial at "x" is faster than using Neville's
      algorithm every time.

   */
   void lagrange( double coef[], double x )
   {
      double x1 = 0;
      double x2 = 1;
      double x3 = 2;
      double x4 = 3;
    
      coef[0] = ((x - x2)* (x - x3) * (x - x4)) / ((x1 - x2) * (x1 - x3) * (x1 - x4));
      coef[1] = ((x - x1)* (x - x3) * (x - x4)) / ((x2 - x1) * (x2 - x3) * (x2 - x4));
      coef[2] = ((x - x1)* (x - x2) * (x - x4)) / ((x3 - x1) * (x3 - x2) * (x3 - x4));
      coef[3] = ((x - x1)* (x - x2) * (x - x3)) / ((x4 - x1) * (x4 - x2) * (x4 - x3));
   } // lagrange

   /**
      Calculate the polynomial interpolation for a point <tt>x</tt>.
     
      The function is passed four data points in <tt>vals</tt> and
      a four value coefficient array, containing the coefficients for
      the polymomial at point <tt>x</tt> (for example 1.5).  This function
      calculates

      <pre>
      f(x) = c<sub>0</sub>y<sub>0</sub> + c<sub>1</sub>y<sub>1</sub> + c<sub>2</sub>y<sub>2</sub> + c<sub>2</sub>y<sub>2</sub>;
      </pre>

   */
   double lagrangeCalc( double vals[], double coef[] )
   {
      double rslt = (vals[0] * coef[0]) + 
         (vals[1] * coef[1]) + 
         (vals[2] * coef[2]) + 
         (vals[3] * coef[3]);

      return rslt;
   } // lagrangeCalc

   /**
      Print the polynomial coefficients.  This is used for debugging.
   */
   void pr_coef( double coef[] ) 
   {
      printf("{%7.4f, %7.4f, %7.4f, %7.4f}\n",
             coef[0], coef[1], coef[2], coef[3] ); 
   }


   public:

   /**
      The interpolation class constructor calculates the coefficients for
      a four point polynomial interpolated at -0.5, 0.5, 1.5, 2.5 and 3.5.
   */
   interpolation() 
   {
      lagrange( coef_minus2,  -2.0 );
      lagrange( coef_minus1,  -1.0 );
      lagrange( coef_minus05, -0.5 );
      lagrange( coef_05,       0.5 );
      lagrange( coef_15,       1.5 );
      lagrange( coef_25,       2.5 );
      lagrange( coef_35,       3.5 );
      lagrange( coef_40,       4.0 );
      lagrange( coef_50,       5.0 );
   } // interpolation class constructor

   /**
      For a point <tt>x</tt> calculate the polynomial interpolation
      from the set of points in <tt>vals</tt>.  The argument <tt>n</tt>
      is the size of the time series (averages).  Note that as
      the wavelet transform calculation recursively proceeds, <tt>n</tt>
      decreases by a factor of two each time (e.g., 32, 16, 8...).
   */
   double pointCalc( double vals[], const double x, const int n )
   {
      double rslt = 0.0;

      if (n == 2) {
         if (x == -0.5) {
            rslt = (vals[0] * (-0.5)) + (vals[1] * 1.5);
         }
         else if (x == 0.5) {
            // interval betwee 0 .. 1
            rslt = (vals[0] * 0.5) + (vals[1] * 0.5);
         }
         else if (x == 1.5) {
            // interval between 1 .. 2
            rslt = (vals[0] * 1.5) + (vals[1] * (-0.5));
         }
         else
            printf("pointCalc, n = 2: can't calculate point at that x = %7.4f\n", x);
      }
      else { // n >= 4
         int ix = (int)(x * 10);
         switch (ix) {
            case -20: // -2.0
            {
               rslt = lagrangeCalc( vals, coef_minus2 );
            }
            break;
            case -10: // -1.0
            {
               rslt = lagrangeCalc( vals, coef_minus1 );
            }
            break;
            case -5: // -0.5
            {
               rslt = lagrangeCalc( vals, coef_minus05 );
            }
            break;
            case 5: // 0.5
            {
               rslt = lagrangeCalc( vals, coef_05 );
            }
            break;
            case 15: // 1.5
            {
               rslt = lagrangeCalc( vals, coef_15 );
            }
            break;
            case 25: // 2.5
            {
               rslt = lagrangeCalc( vals, coef_25 );
            }
            break;
            case 35: // 3.5
            {
               rslt = lagrangeCalc( vals, coef_35 );
            }
            break;
            case 40: // 4.0
            {
               rslt = lagrangeCalc( vals, coef_40 );
            }
            break;
            case 50: // 5.0
            {
               rslt = lagrangeCalc( vals, coef_50 );
            }
            break;
            default:
            {
               printf("pointCalc, n >= 4: can't calculate point at that x = %7.4f\n", x);
            }
            break;
         } // switch
      }

      return rslt;
   } // pointCalc

};  // class interpolation