inplace_haar.java

java的小波分析程序

package wavelets;

import wavelet_util.*;

/**
 *
  
<h4>
   Copyright and Use
</h4>

<p>
   You may use this source code without limitation and without
   fee as long as you include:
</p>
<blockquote>
     This software was written and is copyrighted by Ian Kaplan, Bear
     Products International, www.bearcave.com, 2001.
</blockquote>
<p>
   This software is provided "as is", without any warrenty 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.
<p>
   Please send any bug fixes or suggested source changes to:
<pre>
     [e-mail]
</pre>

<p>
   To generate the documentation for the <tt>wavelets</tt> package
   using Sun Microsystem's <tt>javadoc</tt> use the command
<pre>
        javadoc -private wavelets
</pre>

<p>
   The inplace_haar class calculates an in-place Haar wavelet
   transform.  By in-place it's ment that the result occupies the same
   array as the data set on which the Haar transform is calculated.

<p>
   The Haar wavelet calculation is awkward when the data values are
   not an even power of two.  This is especially true for the in-place
   Haar.  So here we only support data that falls into an even power
   of two.

<p>
   The sad truth about computation is that the time-space tradeoff
   is an iron rule.  The Haar in-place wavelet transform is more 
   memory efficient, but it also takes more computation.

<p>
   The algorithm used here is from <i>Wavelets Made Easy</i>
   by Yves Nievergelt, section 1.4.  The in-place transform
   replaces data values when Haar values and coefficients.
   This algorithm uses a butterfly pattern, where the indices
   are calculated by the following:

<pre>
   for (l = 0; l < log<sub>2</sub>( size ); l++) {
     for (j = 0; j < size; j++) {
        <b>a</b><sub>j</sub> = 2<sup>l</sup> * (2 * j);
        <b>c</b><sub>j</sub> = 2<sup>l</sup> * ((2 * j) + 1);
        if (c<sub>j</sub> >= size)
	  break;
     } <i>// for j</i>
   } <i>// for l</i>
</pre>

<p>
   If there are 16 data elements (indexed 0..15), these loops will
   generate the butterfly index pattern shown below, where the first
   element in a pair is <b>a</b><sub>j</sub>, the Haar value and the
   second element is <b>c</b><sub>j</sub>, the Haar coefficient.

<pre>
{0, 1} {2, 3} {4, 5} {6, 7} {8, 9} {10, 11} {12, 13} {14, 15}
{0, 2} {4, 6} {8, 10} {12, 14}
{0, 4} {8, 12}
{0, 8}
</pre>

<p>
  Each of these index sets represents a Haar wavelet frequency (here
  they are listed from the highest frequency to the lowest).

  @author Ian Kaplan

 */
public class inplace_haar extends wavelet_base {
  /** result of calculating the Haar wavelet */
  private double[] wavefx;
  /** initially false: true means wavefx is ordered by frequency */
  boolean isOrdered = false;

  /**
    Set the wavefx reference variable to the data vector.  Also,
    initialize the isOrdered flag to false.  This indicates that the
    Haar coefficients have not been calculated and ordered by
    frequency.
   */
  public void setWavefx( double[] vector )
  {
    if (vector != null) {
      wavefx = vector;
      isOrdered = false;
    }
  }

  public void setIsOrdered() { isOrdered = true; }

  /**
   *
<p>
   Calculate the in-place Haar wavelet function.  The
   data values are overwritten by the coefficient result, which
   is pointed to by a local reference (<tt>wavefx</tt>).

<p>
   The in-place Haar transform leaves the coefficients in a butterfly
   pattern.  The Haar transform calculates a Haar value
   (<tt><b>a</b><tt>) and a coefficient (<tt>c</tt>) from the forumla
   shown below.

<pre>
  <b>a</b><sub>i</sub> = (v<sub>i</sub> + v<sub>i+1</sub>)/2
   c<sub>i</sub> = (v<sub>i</sub> - v<sub>i+1</sub>)/2
</pre>

<p>

   In the in-place Haar transform the values for <tt><b>a</b></tt> and
   <tt>c</tt> replace the values for v<sub>i</sub> and
   v<sub>i+1</sub>. Subsequent passes calculate new <tt><b>a</b></tt>
   and <tt>c</tt> values from the previous
   <tt><b>a</b><sub>i</sub></tt> and <tt><b>a</b><sub>i+1</sub></tt>
   values.  The produces the butterfly pattern outlined below.

<pre>

v<sub>0</sub> v<sub>1</sub> v<sub>2</sub> v<sub>3</sub> v<sub>4</sub> v<sub>5</sub> v<sub>6</sub> v<sub>7</sub>

<b>a</b><sub>0</sub> c<sub>0</sub> <b>a</b><sub>0</sub> c<sub>0</sub> <b>a</b><sub>0</sub> c<sub>0</sub> <b>a</b><sub>0</sub> c<sub>0</sub>

<b>a</b><sub>1</sub> c<sub>0</sub> c<sub>1</sub> c<sub>0</sub> <b>a</b><sub>1</sub> c<sub>0</sub> c<sub>1</sub> c<sub>0</sub>

<b>a</b><sub>2</sub> c<sub>0</sub> c<sub>1</sub> c<sub>0</sub> c<sub>2</sub> c<sub>0</sub> c<sub>1</sub> c<sub>0</sub>

</pre>

<p> 
   For example, Haar wavelet the calculation with the data set 
   {3, 1, 0, 4, 8, 6, 9, 9} is shown below.  Bold type
   denotes an <b>a</b> value which will be used in the 
   next sweep of the calculation.

<pre>
3  1  0  4  8  6  9  9

<b>2</b>  1  <b>2</b> -2  <b>7</b>  1  <b>9</b>  0

<b>2</b>  1  0 -2  <b>8</b>  1 -1  0

<b>5</b>  1  0 -2 -3  1 -1  0
</pre>


   @param values
      An array of double data values from which the
      Haar wavelet function will be calculated.  The
      number of values must be a power of two.

   */
  public void wavelet_calc( double[] values ) {
    int len = values.length;
    setWavefx( values );

    if (len > 0) {
      byte log = binary.log2( len );
      
      len = binary.power2( log );  // calculation must be on 2 ** n values

      for (byte l = 0; l < log; l++) {
	int p = binary.power2( l );
	for (int j = 0; j < len; j++) {
	  int a_ix = p * (2 * j);
	  int c_ix = p * ((2 * j) + 1);

	  if (c_ix < len) {
	    double a = (values[a_ix] + values[c_ix]) / 2;
	    double c = (values[a_ix] - values[c_ix]) / 2;
	    values[a_ix] = a;
	    values[c_ix] = c;
	  }
	  else {
	    break;
	  }
	} // for j
      } // for l
    }
  } // wavelet_calc


  /**
   *

    Recursively calculate the Haar spectrum, replacing data in the
    original array with the calculated averages.

   */
  private void spectrum_calc(double[] values, 
			     int start, 
			     int end ) {
    int j;
    int newEnd;
    if (start > 0) {
      j = start-1;
      newEnd = end >> 1;
    }
    else {
      j = end-1;
      newEnd = end;
    }

    for (int i = end-1; i > start; i = i - 2, j--) {
      values[j] = (values[i-1] + values[i]) / 2;
    } // for
    

    if (newEnd > 1) {
      int newStart = newEnd >> 1;
      spectrum_calc(values, newStart, newEnd);
    }
  } // spectrum_calc


  /**
   *

<p>
   Calculate the Haar wavelet spectrum

<p> 
   Wavelet calculations sample a signal via a window.  In the case of
   the Haar wavelet, this window is a rectangle.  The signal is
   sampled in passes, using a window of a wider width for each pass.
   Each sampling can be thought of as generating a spectrum of the
   original signal at a particular resolution.

<p>
   In the case of the Haar wavelet, the window is initially two values
   wide.  The first spectrum has half as many values as the original
   signal, where each new value is the average of two values from 
   the original signal.

<p>
   The next spectrum is generated by increasing the window
   width by a factor of two and averaging four elements.
   The third spectrum is generated by increasing the
   window size to eight, etc...  Note that each  of these
   averages can be formed from the previous average.

<p>
   For example, if the initial data set is
<pre>
    { 32.0, 10.0, 20.0, 38.0, 37.0, 28.0, 38.0, 34.0,
      18.0, 24.0, 18.0,  9.0, 23.0, 24.0, 28.0, 34.0 }
</pre>
<p>
    The first spectrum is constructed by averaging
    elements {0,1}, {2,3}, {4,5} ...
</p>
<pre>
    {21, 29, 32.5, 36, 21, 13.5, 23.5, 31} <i>1<sup>st</sup> spectrum</i>
</pre>
<p>
    The second spectrum is constructed by averaging elements
    averaging elements {0,1}, {2,3} in the first spectrum:
</p>
<pre>
    {25, 34.25, 17.25, 27.25}           <i>2<sup>nd</sup> spectrum</i>

    {29.625, 22.25}                     <i>3<sup>ed</sup> spectrum</i>

    {25.9375}                           <i>4<sup>th</sup> spectrum</i>
</pre>
<p>
    Note that we can calculate the Haar spectrum "in-place", by
    replacing the original data values with the spectrum values:
<pre>
    {0, 
     25.9375, 
     29.625, 22.25, 
     25, 34.25, 17.25, 27.25,
     21, 29, 32.5, 36, 21, 13.5, 23.5, 31}
</pre>
<p>
    The spetrum is ordered by increasing frequency.
    This is the same ordering used for the Haar coefficients.
    Keeping to this ordering allows the same code to be applied
    to both the Haar spectrum and a set of Haar coefficients.
</p>
<p>
    This function will destroy the original data.
    When the Haar spectrum is calculated information is lost.
    For example, without the Haar coefficients, which provide the
    difference between the two numbers that form the average, there
    may be several numbers which satisify the equation
<pre>
   <b>a</b><sub>i</sub> = (v<sub>j</sub> + v<sub>j+1</sub>)/2
</pre>

<p>
  For 2<sup>n</sup> initial elements, there will be
  2<sup>n</sup> - 1 results.  For example:
</p>
<pre>
    512 : initial length
    256 : 1st spectrum
    128 : 2nd spectrum
     64 : 3ed spectrum
     32 : 4th spectrum
     16 : 5th spectrum
      8 : 6th spectrum
      4 : 7th spectrum
      2 : 8th spectrum
      1 : 9th spectrum (overall average)
</pre>

<p>
  Since this is an in-place algorithm, the result is
  returned in the values argument.
</p>


   */
  public void wavelet_spectrum( double[] values ) {
    int len = values.length;

    if (len > 0) {
      setWavefx( values );
      byte log = binary.log2( len );
      
      len = binary.power2( log );  // calculation must be on 2 ** n values

      int srcStart = 0;
      spectrum_calc(values, srcStart, len);
      values[0] = 0;
    }
  } // wavelet_vals


  /**
   *

     Print the result of the Haar wavelet function.

   */
  public void pr() {
    if (wavefx != null) {
      int len = wavefx.length;

      System.out.print("{");
      for (int i = 0; i < len; i++) {
	System.out.print( wavefx[i] );
	if (i < len-1)
	  System.out.print(", ");
      }
      System.out.println("}");
    }
  } // pr


  /**
   *
<p>
     Print the Haar value and coefficient showing the 
     ordering.  The Haar value is printed first, followed
     by the coefficients in increasing frequency.  An
     example is shown below.  The Haar value is shown in
     bold.  The coefficients are in normal type.

<p>
    Data set
<pre>
    { 32.0, 10.0, 20.0, 38.0,
      37.0, 28.0, 38.0, 34.0,
      18.0, 24.0, 18.0,  9.0, 
      23.0, 24.0, 28.0, 34.0 }
</pre>
<p>
    Ordered Haar transform:
<pre>
   <b>25.9375</b>
   3.6875
   -4.625 -5.0
   -4.0 -1.75 3.75 -3.75