lifting.h

小波提升格式的源代码

   where the ordering is defined by the sub-classes.

 */
class liftingScheme {
   private:
    interpolation fourPtInterp;
   
   /** Enumerator for add and subtract operators */
   typedef enum { bad_dir, forward, inverse } direction;


  /**
    <p>
    Copy four points or <i>N</i> (which ever is less) data points from
    <i>vec</i> into <i>d</i> These points are the "known" points used
    in the polynomial interpolation.
    </p>

    @param vec[] the input data set on which the wavelet is calculated
    @param d[] an array into which <i>N</i> data points, starting at
               <i>start</i> are copied.
    @param N the number of polynomial interpolation points
    @param start the index in <i>vec</i> from which copying starts
   */
  void fill( DoubleVec& vec, double d[], int N, int start )
  {
    int n = 4;
    if (n > N)
      n = N;
    int end = start + n;
    int j = 0;

    for (int i = start; i < end; i++) {
      d[j] = vec[i];
      j++;
    }
  } // fill


  /**
    <p>
    Predict an odd point from the even points, using 4-point
    polynomial interpolation.
    </p>
    <p>
    The four points used in the polynomial interpolation are
    the even points.  We pretend that these four points
    are located at the x-coordinates 0,1,2,3.  The first 
    odd point interpolated will be located between the first
    and second even point, at 0.5.  The next N-3 points are
    located at 1.5 (in the middle of the four points).
    The last two points are located at 2.5 and 3.5.
    For complete documentation see
    </p>
    <pre>
  <a href="http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html">
  http://www.bearcave.com/misl/misl_tech/wavelets/lifting/index.html</a>
    </pre>

    <p>
    The difference between the predicted (interpolated) value
    and the actual odd value replaces the odd value in the
    forward transform.
    </p>

    <p>
    As the recursive steps proceed, N will eventually be 4 and
    then 2.  When N = 4, linear interpolation is used.
    When N = 2, Haar interpolation is used (the prediction for
    the odd value is that it is equal to the even value).
    </p>

    @param vec the input data on which the forward or inverse
               transform is calculated.
    @param N the area of vec over which the transform is calculated
    @param direction forward or inverse transform

   */
  void interp( DoubleVec& vec, const int N, const direction dir )
  {
    int half = N >> 1;
    double d[4];

    for (int i = 0; i < half; i++) {
      double predictVal;

      if (i == 0) {
	if (half == 1) {
	  // e.g., N == 2, and we use Haar interpolation
	  predictVal = vec[0];
	}
	else {
	  fill( vec, d, N, 0 );
	  predictVal = fourPtInterp.interpPoint( 0.5, half, d );
	}
      }
      else if (i == 1) {
	predictVal = fourPtInterp.interpPoint( 1.5, half, d );
      }
      else if (i == half-2) {
	predictVal = fourPtInterp.interpPoint( 2.5, half, d );
      }
      else if (i == half-1) {
	predictVal = fourPtInterp.interpPoint( 3.5, half, d );
      }
      else {
	fill( vec, d, N, i-1);
	predictVal = fourPtInterp.interpPoint( 1.5, half, d );
      }

      int j = i + half;
      if (dir == forward) {
	vec[j] = vec[j] - predictVal;
      }
      else if (dir == inverse) {
	vec[j] = vec[j] + predictVal;
      }
      else {
	printf("interp: bad direction value\n");
	break;
      }
    } // for
  } // interp


   /**
      For the first <tt>N</tt> elements of <tt>ts</tt>, move the odd elements
      to the second half of the vector, leaving the even elements in the first
      half.

      In terms of the wavelet transform, this divides the vector
      into an average region and a coefficient (difference) region.

      For example, if we label the even elements with <i>a</i> and the 
      odd elements with <i>b</i> the input vector will contain

<pre>
  a<sub>0</sub> b<sub>0</sub> a<sub>1</sub> b<sub>1</sub> a<sub>2</sub> b<sub>2</sub> a<sub>3</sub> b<sub>3</sub> 
</pre>

      The split function will move the odd elements into the second half
      of the array, resulting in

<pre>
  a<sub>0</sub> a<sub>2</sub> a<sub>3</sub> a<sub>4</sub> b<sub>0</sub> b<sub>2</sub> b<sub>3</sub> b<sub>4</sub> 
</pre>

   */
   void split( DoubleVec& ts, const int N )
   {
      int start = 1;
      int end = N-1;

      while (start < end) {
         int i;
         double tmp;
         for (i = start; i < end; i = i + 2) {
            tmp = ts[i];
            ts[i] = ts[i+1];
            ts[i+1] = tmp;
         }
         start = start + 1;
         end = end - 1;
      }
   } // split


   /**
      Merge the odd and even values.  The is the inverse of the
      split.

      If the even values are labeled with <i>a</i> and the odd values
      are labeled with <i>b</i>, the input vector will contain

<pre>
  a<sub>0</sub> a<sub>2</sub> a<sub>3</sub> a<sub>4</sub> b<sub>0</sub> b<sub>2</sub> b<sub>3</sub> b<sub>4</sub> 
</pre>

      The merge function will re-order the odd and even element in the
      vector:

<pre>
  a<sub>0</sub> b<sub>0</sub> a<sub>1</sub> b<sub>1</sub> a<sub>2</sub> b<sub>2</sub> a<sub>3</sub> b<sub>3</sub> 
</pre>
    */
   void merge( DoubleVec& ts, int N)
   {
      int half = N >> 1;
      int start = half-1;
      int end = half;
      int i;
      double tmp;

      while (start > 0) {
         for (i = start; i < end; i = i + 2) {
            tmp = ts[i];
            ts[i] = ts[i+1];
            ts[i+1] = tmp;
         }
         start = start - 1;
         end = end + 1;
      }
   } // merge


   /**
      Predict step of the wavelet Lifting Scheme.

      The term "predict" comes from <i>Building Your Own Wavelets
      at Home</i>.

      <b>transform</b>

      In the wavelet transform, calculate the difference between even
      and odd element pairs.  This is a slightly modified version of
      the classic haar difference.  If the even elements are labeled
      as <i>a</i> and the odd elements are labeled as <i>b</i> (where
      we start counting at zero) the difference is simply

<pre>
       d<sub>i</sub> = b<sub>i</sub> - a<sub>i</sub>
</pre>

      Since an "in-place" algorithm is used, where we update the
      odd elements with the difference, we get
<pre>
       b<sub>i</sub> = b<sub>i</sub> - a<sub>i</sub>
</pre>

      <b>inverse transform</b>

      Reverse the transform predict step by adding the average
      (even element) to the difference (odd element).

    */
   void predict( DoubleVec& ts, const int N, const direction dir)
   {
      int i;
      int half = N >> 1;
      int j = 0;
      for (i = half; i < N; i++) {
         if (dir == forward) { // forward transform stage
	   ts[i] = ts[i] - ts[j];
         }
         else if (dir == inverse ) { // inverse transform stage
	   ts[i] = ts[i] + ts[j];
         }
	 else {
	   printf("predict: bad direction\n");
	   break;
	 }
         j++;
      }
   } // predict


   /**
      The update step of the wavelet Lifting Scheme

      <b>transform</b>

      In the Lifting Scheme transform the update step follows
      the predict step.  After the predict step, the differences
      have been calculated in-place writing over the even (b)
      values.  The update step calculates the Haar average using
      an even/odd pair.  The average is placed in the even member
      of the pair.

    */
   void update( DoubleVec& ts, int N, const direction dir)
   {
      int i;
      int half = N >> 1;
      int j = half;
      for (i = 0; i < half; i++) {
         if (dir == forward) {  // forward transform stage
            ts[i] = ts[i] + (ts[j]/2.0);
         }
         else if (dir == inverse) { // inverse transform step
            ts[i] = ts[i] - (ts[j]/2.0);
         }
	 else {
	   printf("update: bad direction value\n");
	   break;
	 }
         j++;
      } // for
   } // update


 public:
   /**
      Wavelet lifting scheme transform.

      The wavelet lifting scheme is an in-place wavelet algorithm.
      A time series of N elements is passed to the <tt>transform</tt>
      function.  The result of the wavelet lifting scheme is calculated
      in place (without any array temporaries) in the <tt>ts</tt>
      DoubleVec.

    */
   void transform( DoubleVec& ts )
   {
     const int N = ts.size();
     int n;
     for (n = N; n > 1; n = n >> 1) {
       split( ts, n );
       predict( ts, n, forward );
       update( ts, n, forward );
       interp( ts, n, forward );
     } // for
   } // transform

   /**
      Wavelet lifting Scheme inverse transform.

      Like the Lifting Scheme transform, the inverse transform is an
      in-place algorithm.  The result of the transform is passed to
      the inverse transform, which calculates the time series from
      the time series average and the coefficients.

    */
   void invTransform( DoubleVec& ts )
   {
     const int N = ts.size();
     int n;
     for (n = 2; n <= N; n = n << 1) {
       interp( ts, n, inverse );
       update( ts, n, inverse );
       predict( ts, n, inverse );
       merge( ts, n );
     }
   } // invTransform

}; // class liftingScheme