inplace_haar.java

java的小波分析程序

   11.0 -9.0 4.5 2.0 -3.0 4.5 -0.5 -3.0
</pre>

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

      if (len > 0) {
	System.out.println(wavefx[0]);

	int num_in_freq = 1;
	int cnt = 0;
	for (int i = 1; i < len; i++) {
	  System.out.print(wavefx[i] + " ");
	  cnt++;
	  if (cnt == num_in_freq) {
	    System.out.println();
	    cnt = 0;
	    num_in_freq = num_in_freq * 2;
	  }
	}
      }
    }
  } // pr_ordered


  /**
   *
<p>

     Order the result of the in-place Haar wavelet function,
     referenced by wavefx.  As noted above in the documentation for
     <tt>wavelet_calc()</tt>, the in-place Haar transform leaves the
     coefficients in a butterfly pattern.  This can be awkward for
     some calculations.  The <tt>order</tt> function orders the
     coefficients by frequency, from the lowest frequency to the
     highest.  The number of coefficients for each frequency band
     follow powers of two (e.g., 1, 2, 4, 8 ... 2<sup>n</sup>).
     An example of the coefficient sort performed by the <tt>order()</tt>
     function is shown below;

<pre>
before: <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>

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

</pre>

<p>
     The results in the same ordering as the ordered Haar wavelet
     transform.

<p>
     If the number of elements is 2<sup>n</sup>, then the largest
     number of coefficients will be 2<sup>n-1</sup>.  Each of the
     coefficients in the largest group is separated by one element
     (which contains other coefficients).  This algorithm pushes these
     together so they are not separated.  These coefficients now make
     up half of the array.  The remaining coefficients take up the
     other half.  The next frequency down is also separated by one
     element.  These are pushed together taking up half of the half.
     The algorithm keeps going until only one coefficient is left.

<p>
     As with wavelet_calc above, this algorithm assumes that
     the number of values is a power of two.

   */
  public void order() {
    if (wavefx != null) {

      int half = 0;
      for (int len = wavefx.length; len > 1; len = half) {
	half = len / 2;
	int skip = 1;
	for (int dest = len - 2; dest >= half; dest--) {
	  int src = dest - skip;
	  double tmp = wavefx[src];
	  for (int i = src; i < dest; i++)
	    wavefx[i] = wavefx[i+1];
	  wavefx[dest] = tmp;
	  skip++;
	} // for dest
      } // for len

      isOrdered = true;
    } // if
  } // order()


  /**
   *
<p>
     Regenerate the data from the Haar wavelet function.

<p>

     There is no information loss in the Haar function.  The original
     data can be regenerated by reversing the calculation.  Given
     a Haar value, <b>a</b> and a coefficient <b>c</b>, two Haar
     values can be generated

<pre>
        <b>a</b><sub>i</sub>  = a + c;
        <b>a</b><sub>i+1</sub> = a - c;
</pre>

<p>
     The transform is calculated from the low frequency coefficients
     to the high frequency coefficients.  An example is shown below
     for the result of the ordered Haar transform.  Note that the
     values are in bold and the coefficients are in normal type.

<p>
To regenerate {<a>a</b><sub>1</sub>, <a>a</b><sub>2</sub>, <a>a</b><sub>3</sub>, <a>a</b><sub>4</sub>, <a>a</b><sub>5</sub>, <a>a</b><sub>6</sub>, <a>a</b><sub>7</sub>, <a>a</b><sub>8</sub>} from
<pre>
<b>a</b><sub>1</sub>
c<sub>1</sub>
c<sub>2</sub> c<sub>3</sub>
c<sub>4</sub> c<sub>5</sub> c<sub>6</sub> c<sub>7</sub>
</pre>
<p>
   The inverse Haar transform is applied:

<pre>
<b>a</b><sub>1</sub> = <b>a</b><sub>1</sub> + c<sub>1</sub>
<b>a</b><sub>2</sub> = <b>a</b><sub>1</sub> - c<sub>1</sub>

<b>a</b><sub>1</sub> = <b>a</b><sub>1</sub> + c<sub>2</sub>
<b>a</b><sub>2</sub> = <b>a</b><sub>1</sub> - c<sub>2</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> + c<sub>3</sub>
<b>a</b><sub>4</sub> = <b>a</b><sub>2</sub> - c<sub>3</sub>

<b>a</b><sub>1</sub> = <b>a</b><sub>1</sub> + c<sub>4</sub>
<b>a</b><sub>2</sub> = <b>a</b><sub>1</sub> - c<sub>4</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> + c<sub>5</sub>
<b>a</b><sub>4</sub> = <b>a</b><sub>2</sub> - c<sub>5</sub>
<b>a</b><sub>5</sub> = <b>a</b><sub>3</sub> + c<sub>6</sub>
<b>a</b><sub>6</sub> = <b>a</b><sub>3</sub> - c<sub>6</sub>
<b>a</b><sub>7</sub> = <b>a</b><sub>4</sub> + c<sub>7</sub>
<b>a</b><sub>8</sub> = <b>a</b><sub>4</sub> - c<sub>7</sub>
</pre>
<p>
    For example:

<pre>
<b>5.0</b>
-3.0
0.0 -1.0
1.0 -2.0 1.0 0.0

<b>5.0</b>+(-3), <b>5.0</b>-(-3) = {<b>2</b> <b>8</b>}

<b>2</b>+0, <b>2</b>-0, <b>8</b>+(-1), <b>8</b>-(-1) = {<b>2, 2, 7, 9</b>}

<b>2</b>+1, <b>2</b>-1, <b>2</b>+(-2), <b>2</b>-(-2), <b>7</b>+1, <b>7</b>-1, <b>9</b>+0, <b>9</b>-0 = {3,1,0,4,8,6,9,9}

</pre>

<p>
   By using the butterfly indices the inverse transform can
   also be applied to an unordered in-place haar function.

<p>
   This function checks the to see whether the wavefx array is
   ordered.  If wavefx is ordered the inverse transform described
   above is applied.  If the data remains in the in-order configuration
   an inverse butterfly is applied.  Note that once the inverse
   Haar is calculated the Haar function data will be replaced by
   the original data.

   */
  public void inverse() {
    if (wavefx != null) {
      if (isOrdered) {
	inverse_ordered();
	// Since the signal has been rebuilt from the
	// ordered coefficients, set isOrdered to false
	isOrdered = false;
      }
      else {
	inverse_butterfly();
      }
    }
  } // inverse


  /**
   *
<p>
     Calculate the inverse Haar transform on an ordered
     set of coefficients.
<p>
     See comment above for the <tt>inverse()</tt> method
     for details.
<p>
     The algorithm above uses an implicit temporary.  The
     in-place algorithm is a bit more complex.  As noted
     above, the Haar value and coefficient are replaced
     with the newly calculated values:

<pre>
     t<sub>1</sub> = <b>a</b><sub>1</sub> + c<sub>1</sub>;
     t<sub>2</sub> = <b>a</b><sub>1</sub> - c<sub>1</sub>;
     <b>a</b><sub>1</sub> = t<sub>1</sub>;
     c<sub>1</sub> = t<sub>2</sub>
</pre>

<p>
    As the calculation proceeds and the coefficients are replaced by
    the newly calculated Haar values, the values are out of order.
    This is shown in the below (use <tt>javadoc -private
    wavelets</tt>).  Here each element is numbered with a subscript as
    it should be ordered.  A sort operation reorders these values
    as the calculation proceeds.  The variable <tt>L</tt> is the power of
    two.

<pre>
start: {<b>5.0</b>, -3.0, 0.0, -1.0, 1.0, -2.0, 1.0, 0.0}
[0, 1]
L = 1
{<b>2.0</b><sub>0</sub>, <b>8.0</b><sub>1</sub>, 0.0, -1.0, 1.0, -2.0, 1.0, 0.0}
 sort:

 start: {<b>2.0</b><sub>0</sub>, <b>8.0</b><sub>1</sub>, 0.0, -1.0, 1.0, -2.0, 1.0, 0.0}
[0, 2], [1, 3]
L = 2
{<b>2.0</b><sub>0</sub>, <b>7.0</b><sub>2</sub>, <b>2.0</b><sub>1</sub>, <b>9.0</b><sub>3</sub>, 1.0, -2.0, 1.0, 0.0}
 sort:
exchange [1, 2]
{<b>2.0</b><sub>0</sub>, <b>2.0</b><sub>1</sub>, <b>7.0</b><sub>2</sub>, <b>9.0</b><sub>3</sub>, 1.0, -2.0, 1.0, 0.0}

 start: {2.0, 2.0, 7.0, 9.0, 1.0, -2.0, 1.0, 0.0}
[0, 4], [1, 5], [2, 6], [3, 7]
L = 4
{<b>3.0</b><sub>0</sub>, <b>0.0</b><sub>2</sub>, <b>8.0</b><sub>4</sub>, <b>9.0</b><sub>6</sub>, <b>1.0</b><sub>1</sub>, <b>4.0</b><sub>3</sub>, <b>6.0</b><sub>5</sub>, <b>9.0</b><sub>7</sub></b>}
 sort:
exchange [1, 4]
{<b>3.0</b><sub>0</sub>, <b>1.0</b><sub>1</sub>, <b>8.0</b><sub>4</sub>, <b>9.0</b><sub>6</sub>, <b>0.0</b><sub>2</sub>, <b>4.0</b><sub>3</sub>, <b>6.0</b><sub>5</sub>, <b>9.0</b><sub>7</sub>}
exchange [2, 5]
{<b>3.0</b><sub>0</sub>, <b>1.0</b><sub>1</sub>, <b>4.0</b><sub>3</sub>, <b>9.0</b><sub>6</sub>, <b>0.0</b><sub>2</sub>, <b>8.0</b><sub>4</sub>, <b>6.0</b><sub>5</sub>, <b>9.0</b><sub>7</sub>}
exchange [3, 6]
{<b>3.0</b><sub>0</sub>, <b>1.0</b><sub>1</sub>, <b>4.0</b><sub>3</sub>, <b>6.0</b><sub>5</sub>, <b>0.0</b><sub>2</sub>, <b>8.0</b><sub>4</sub>, <b>9.0</b><sub>6</sub>, <b>9.0</b><sub>7</sub>}
exchange [2, 4]
{<b>3.0</b><sub>0</sub>, <b>1.0</b><sub>1</sub>, <b>0.0</b><sub>2</sub>, <b>6.0</b><sub>5</sub>, <b>4.0</b><sub>3</sub>, <b>8.0</b><sub>4</sub>, <b>9.0</b><sub>6</sub>, <b>9.0</b><sub>7</sub>}
exchange [3, 5]
{<b>3.0</b><sub>0</sub>, <b>1.0</b><sub>1</sub>, <b>0.0</b><sub>2</sub>, <b>8.0</b><sub>4</sub>, <b>4.0</b><sub>3</sub>, <b>6.0</b><sub>5</sub>, <b>9.0</b><sub>6</sub>, <b>9.0</b><sub>7</sub>}
exchange [3, 4]
{<b>3.0</b><sub>0</sub>, <b>1.0</b><sub>1</sub>, <b>0.0</b><sub>2</sub>, <b>4.0</b><sub>3</sub>, <b>8.0</b><sub>4</sub>, <b>6.0</b><sub>5</sub>, <b>9.0</b><sub>6</sub>, <b>9.0</b><sub>7</sub>}
</pre>   
 
  ********/

  private void inverse_ordered() {
    int len = wavefx.length;
    
    for (int L = 1; L < len; L = L * 2) {
      
      int i;
      
      // calculate
      for (i = 0; i < L; i++) {
	int a_ix = i;
	int c_ix = i + L;
	double a1        = wavefx[a_ix] + wavefx[c_ix];
	double a1_plus_1 = wavefx[a_ix] - wavefx[c_ix];
	wavefx[a_ix] = a1;
	wavefx[c_ix] =a1_plus_1;
      } // for i
      
      // sort
      int offset = L-1;
      for (i = 1; i < L; i++) {
	for (int j = i; j < L; j++) {
	  double tmp = wavefx[j];
	  wavefx[j] = wavefx[j+offset];
	  wavefx[j+offset] = tmp;
	} // for j
	offset = offset - 1;
      } // for i
      
    } // for L
  } // inverse_ordered


  /**
   *
<p>
    Calculate the inverse Haar transform on the result of
    the in-place Haar transform.

<p>
    The inverse butterfly exactly reverses in-place Haar
    transform.  Instead of generating coefficient values
    (<tt>c<sub>i</sub><tt>), the inverse butterfly calculates
    Haar values (<tt><b>a</b><sub>i</sub>) using the 
    equations:
<pre>
        <b>new_a</b><sub>i</sub> = <b>a</b><sub>i</sub> + c<sub>i</sub>
        <b>new_a</b><sub>i+1</sub> = <b>a</b><sub>i</sub> - c<sub>i</sub>
</pre>

<pre>
<b>a</b><sub>0</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>

<b>a</b><sub>1</sub> = <b>a</b><sub>0</sub> + c<sub>2</sub>
<b>a</b><sub>1</sub> = <b>a</b><sub>0</sub> - c<sub>2</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> = <b>a</b><sub>1</sub> + c<sub>1</sub>
<b>a</b><sub>2</sub> = <b>a</b><sub>1</sub> - c<sub>1</sub>
<b>a</b><sub>2</sub> = <b>a</b><sub>1</sub> + c<sub>1</sub>
<b>a</b><sub>2</sub> = <b>a</b><sub>1</sub> - c<sub>1</sub>

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

<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> + c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> - c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> + c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> - c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> + c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> - c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> + c<sub>0</sub>
<b>a</b><sub>3</sub> = <b>a</b><sub>2</sub> - c<sub>0</sub>

<b>a</b><sub>3</sub> <b>a</b><sub>3</sub> <b>a</b><sub>3</sub> <b>a</b><sub>3</sub> <b>a</b><sub>3</sub> <b>a</b><sub>3</sub> <b>a</b><sub>3</sub> <b>a</b><sub>3</sub> 

</pre>

<p>
A numeric example is shown below.

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

<b>a</b><sub>1</sub> = 5 + (-3) = 2
<b>a</b><sub>1</sub> = 5 - (-3) = 8

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

<b>a</b><sub>2</sub> = 2 + 0    = 2
<b>a</b><sub>2</sub> = 2 - 0    = 2
<b>a</b><sub>2</sub> = 8 + (-1) = 7
<b>a</b><sub>2</sub> = 8 - (-1) = 9

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

<b>a</b><sub>3</sub> = 2 + 1    = 3
<b>a</b><sub>3</sub> = 2 - 1    = 1
<b>a</b><sub>3</sub> = 2 + (-2) = 0
<b>a</b><sub>3</sub> = 2 - (-2) = 4
<b>a</b><sub>3</sub> = 7 + 1    = 8
<b>a</b><sub>3</sub> = 7 - 1    = 6
<b>a</b><sub>3</sub> = 9 + 0    = 9
<b>a</b><sub>3</sub> = 9 - 0    = 9

<b>3</b><sub>3</sub>  <b>1</b><sub>3</sub>  <b>0</b><sub>3</sub>  <b>4</b><sub>3</sub>  <b>8</b><sub>3</sub>  <b>6</b><sub>3</sub>  <b>9</b><sub>3</sub>  <b>9</b><sub>3</sub>

</pre>

<p>
   Note that the inverse_butterfly function is faster than
   the inverse_ordered function, since data does not have
   to be reshuffled during the calculation.

   */
  private void inverse_butterfly() {
    int len = wavefx.length;
    
    if (len > 0) {
      byte log = binary.log2( len );
      
      len = binary.power2( log );  // calculation must be on 2 ** n values
      
      for (byte l = (byte)(log-1); l >= 0; 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 = wavefx[a_ix] + wavefx[c_ix];
	    double c = wavefx[a_ix] - wavefx[c_ix];
	    wavefx[a_ix] = a;
	    wavefx[c_ix] = c;
	  }
	  else {
	    break;
	  }
	} // for j
      } // for l
    }
  } // inverse_butterfly


} // inplace_haar