line.h

小波包分解去噪c++源程序

#ifndef _LINE_H_
#define _LINE_H_

#include "liftbase.h"


/**
  Line (with slope) wavelet

  The wavelet Lifting Scheme "line" wavelet approximates the
  data set using a line with with slope (in contrast to the
  Haar wavelet where a line has zero slope is used to approximate
  the data).

  The predict stage of the line wavelet "predicts" that an odd point
  will lie midway between its two neighboring even points.  That is,
  that the odd point will lie on a line between the two adjacent even
  points.  The difference between this "prediction" and the actual odd
  value replaces the odd element.

  The update stage calculates the average of the odd and even 
  element pairs, although the method is indirect, since the 
  predict phase has over written the odd value.

<h4>
   Copyright and Use
</h4>

   You may use this source code without limitation and without
   fee as long as you include:

<blockquote>
     This software was written and is copyrighted by Ian Kaplan, Bear
     Products International, www.bearcave.com, 2001.
</blockquote>

   This software is provided "as is", without any warranty or
   claim as to its usefulness.  Anyone who uses this source code
   uses it at their own risk.  Nor is any support provided by
   Ian Kaplan and Bear Products International.

   Please send any bug fixes or suggested source changes to:

<pre>
     [e-mail]
</pre>

  \author Ian Kaplan

 */
template <class T>
class line : public liftbase<T, double> {

private:

   /**
     Calculate an extra "even" value for the line wavelet algorithm 
     at the end of the data series.  Here we pretend that the
     last two values in the data series are at the x-axis
     coordinates 0 and 1, respectively.  We then need to calculate
     the y-axis value at the x-axis coordinate 2.  This point lies
     on a line running through the points at 0 and 1.

     Given two points, x<sub>1</sub>, y<sub>1</sub> and 
     x<sub>2</sub>, y<sub>2</sub>, where

     <pre>
        x<sub>1</sub> = 0
        x<sub>2</sub> = 1
     </pre>

     calculate the point on the line at x<sub>3</sub>, y<sub>3</sub>,
     where 

     <pre>
        x<sub>3</sub> = 2
     </pre>

     The "two-point equation" for a line given x<sub>1</sub>, 
     y<sub>1</sub> and x<sub>2</sub>, y<sub>2</sub> is

     <pre>
     .          y<sub>2</sub> - y<sub>1</sub>
     (y - y<sub>1</sub>) = -------- (x - x<sub>1</sub>)
     .          x<sub>2</sub> - x<sub>1</sub>
     </pre>

     Solving for y

     <pre>
     .    y<sub>2</sub> - y<sub>1</sub>
     y = -------- (x - x<sub>1</sub>) + y<sub>1</sub>
     .    x<sub>2</sub> - x<sub>1</sub>
     </pre>

     Since x<sub>1</sub> = 0 and x<sub>2</sub> = 1

     <pre>
     .    y<sub>2</sub> - y<sub>1</sub>
     y = -------- (x - 0) + y<sub>1</sub>
     .    1 - 0
     </pre>

     or

     <pre>
     y = (y<sub>2</sub> - y<sub>1</sub>)*x + y<sub>1</sub>
     </pre>

     We're calculating the value at x<sub>3</sub> = 2, so

     <pre>
     y = 2*y<sub>2</sub> - 2*y<sub>1</sub> + y<sub>1</sub>
     </pre>

     or

     <pre>
     y = 2*y<sub>2</sub> - y<sub>1</sub>
     </pre>

    */
   double new_y( double y1, double y2)
   {
     double y = 2 * y2 - y1;
     return y;
   }


protected:

  /**
    Predict phase of line Lifting Scheme wavelet
    
    The predict step attempts to "predict" the value of an
    odd element from the even elements.  The difference
    between the prediction and the actual element is stored
    as a wavelet coefficient.

    The "predict" step takes place after the split step.  The split
    step will move the odd elements (b<sub>j</sub>) to the second half
    of the array, leaving the even elements (a<sub>i</sub>) in the
    first half

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

    The predict step of the line wavelet "predicts" that the
    odd element will be on a line between two even elements.

    <pre>
    b<sub>j+1,i</sub> = b<sub>j,i</sub> - (a<sub>j,i</sub> + a<sub>j,i+1</sub>)/2
    </pre>

    Note that when we get to the end of the data series the odd
    element is the last element in the data series (remember,
    wavelet algorithms work on data series with 2<sup>n</sup>
    elements).  Here we "predict" that the odd element will be
    on a line that runs through the last two even elements.
    This can be calculated by assuming that the last two
    even elements are located at x-axis coordinates 0 and 1,
    respectively.  The odd element will be at 2.  The <i>new_y()</i>
    function is called to do this simple calculation.

   */
  void predict( T& vec, int N, transDirection direction )
  {
    int half = N >> 1;
    double predictVal;

    for (int i = 0; i < half; i++) {
      int j = i + half;
      if (i < half-1) {
	predictVal = (vec[i] + vec[i+1])/2;
      }
      else if (N == 2) {
	predictVal = vec[0];
      }
      else {
	// calculate the last "odd" prediction
	double n_plus1 = new_y( vec[i-1], vec[i] );
	predictVal = (vec[i] + n_plus1)/2;
      }

      if (direction == forward) {
	vec[j] = vec[j] - predictVal;
      }
      else if (direction == inverse) {
	vec[j] = vec[j] + predictVal;
      }
      else {
	printf("predictline::predict: bad direction value\n");
      }
    }
  } // predict

  
  /**
    Update step of the linear interpolation wavelet

    The predict phase works on the odd elements in the second
    half of the array.  The update phase works on the even
    elements in the first half of the array.  The update
    phase attempts to preserve the average.  After the update
    phase is completed the average of the even elements should
    be approximately the same as the average of the input data
    set from the previous iteration.  The result of the update
    phase becomes the input for the next iteration.

    In a Haar wavelet the average that replaces the even element is
    calculated as the average of the even element and its associated
    odd element (e.g., its odd neighbor before the split).  This is
    not possible in the line wavelet since the odd element has been
    replaced by the difference between the odd element and the
    mid-point of its two even neighbors.  As a result, the odd element
    cannot be recovered.

    The value that is added to the even element to preserve the
    average is calculated by the equation shown below.  This
    equation is given in Wim Sweldens' journal articles and his tutorial
    (<i>Building Your Own Wavelets at Home</i>) and in 
    <i>Ripples in Mathematics</i>.  A somewhat more complete derivation
    of this equation is provided in <i>Ripples in Mathematics</i>
    by A. Jensen and A. la Cour-Harbo, Springer, 2001.

    The equation used to calculate the average is shown below for a
    given iteratin <i>i</i>.  Note that the predict phase has already
    completed, so the odd values belong to iteration <i>i+1</i>.

    <pre>
  even<sub>i+1,j</sub> = even<sub>i,j</sub> op (odd<sub>i+1,k-1</sub> + odd<sub>i+1,k</sub>)/4
    </pre>

    There is an edge problem here, when i = 0 and k = N/2 (e.g., there
    is no k-1 element).  We assume that the odd<sub>i+1,k-1</sub> is
    the same as odd<sub>k</sub>.  So for the first element this
    becomes

    <pre>
      (2 * odd<sub>k</sub>)/4
    </pre>

    or

    <pre>
      odd<sub>k</sub>/2
    </pre>

   */
  void update( T& vec, int N, transDirection direction )
  {
    int half = N >> 1;

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

      if (i == 0) {
	val = vec[j]/2.0;
      }
      else {
	val = (vec[j-1] + vec[j])/4.0;
      }
      if (direction == forward) {
	vec[i] = vec[i] + val;
      }
      else if (direction == inverse) {
	vec[i] = vec[i] - val;
      }
      else {
	printf("update: bad direction value\n");
      }
    } // for    
  }

}; // line

#endif