View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.transform;
18  
19  import org.apache.commons.math.FunctionEvaluationException;
20  import org.apache.commons.math.MathRuntimeException;
21  import org.apache.commons.math.analysis.UnivariateRealFunction;
22  
23  /**
24   * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
25   * Transformation of an input vector x to the output vector y.
26   * <p>In addition to transformation of real vectors, the Hadamard transform can
27   * transform integer vectors into integer vectors. However, this integer transform
28   * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
29   * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
30   * vector (1/2, -1/2, 0, 0).</p>
31   * @version $Revision: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 2009) $
32   * @since 2.0
33   */
34  public class FastHadamardTransformer implements RealTransformer {
35  
36      /** {@inheritDoc} */
37      public double[] transform(double f[])
38          throws IllegalArgumentException {
39          return fht(f);
40      }
41  
42      /** {@inheritDoc} */
43      public double[] transform(UnivariateRealFunction f,
44                                double min, double max, int n)
45          throws FunctionEvaluationException, IllegalArgumentException {
46          return fht(FastFourierTransformer.sample(f, min, max, n));
47      }
48  
49      /** {@inheritDoc} */
50      public double[] inversetransform(double f[])
51      throws IllegalArgumentException {
52          return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
53     }
54  
55      /** {@inheritDoc} */
56      public double[] inversetransform(UnivariateRealFunction f,
57                                       double min, double max, int n)
58          throws FunctionEvaluationException, IllegalArgumentException {
59          final double[] unscaled =
60              fht(FastFourierTransformer.sample(f, min, max, n));
61          return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
62      }
63  
64      /**
65       * Transform the given real data set.
66       * <p>The integer transform cannot be inverted directly, due to a scaling
67       * factor it may lead to double results.</p>
68       * @param f the integer data array to be transformed (signal)
69       * @return the integer transformed array (spectrum)
70       * @throws IllegalArgumentException if any parameters are invalid
71       */
72      public int[] transform(int f[])
73          throws IllegalArgumentException {
74          return fht(f);
75      }
76  
77      /**
78       * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
79       * <br>
80       * Requires <b>Nlog2N = n2</b><sup>n</sup> additions.
81       * <br>
82       * <br>
83       * <b><u>Short Table of manual calculation for N=8:</u></b>
84       * <ol>
85       * <li><b>x</b> is the input vector we want to transform</li>
86       * <li><b>y</b> is the output vector which is our desired result</li>
87       * <li>a and b are just helper rows</li>
88       * </ol>
89       * <pre>
90       * <code>
91       * +----+----------+---------+----------+
92       * | <b>x</b>  |    <b>a</b>     |    <b>b</b>    |    <b>y</b>     |
93       * +----+----------+---------+----------+
94       * | x<sub>0</sub> | a<sub>0</sub>=x<sub>0</sub>+x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>+a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>+b<sub>1</sub> |
95       * +----+----------+---------+----------+
96       * | x<sub>1</sub> | a<sub>1</sub>=x<sub>2</sub>+x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>+a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>+b<sub>3</sub> |
97       * +----+----------+---------+----------+
98       * | x<sub>2</sub> | a<sub>2</sub>=x<sub>4</sub>+x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>+a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>+b<sub>5</sub> |
99       * +----+----------+---------+----------+
100      * | x<sub>3</sub> | a<sub>3</sub>=x<sub>6</sub>+x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>+a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>+b<sub>7</sub> |
101      * +----+----------+---------+----------+
102      * | x<sub>4</sub> | a<sub>0</sub>=x<sub>0</sub>-x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>-a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>-b<sub>1</sub> |
103      * +----+----------+---------+----------+
104      * | x<sub>5</sub> | a<sub>1</sub>=x<sub>2</sub>-x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>-a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>-b<sub>3</sub> |
105      * +----+----------+---------+----------+
106      * | x<sub>6</sub> | a<sub>2</sub>=x<sub>4</sub>-x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>-a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>-b<sub>5</sub> |
107      * +----+----------+---------+----------+
108      * | x<sub>7</sub> | a<sub>3</sub>=x<sub>6</sub>-x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>-a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>-b<sub>7</sub> |
109      * +----+----------+---------+----------+
110      * </code>
111      * </pre>
112      * 
113      * <b><u>How it works</u></b>
114      * <ol>
115      * <li>Construct a matrix with N rows and n+1 columns<br>   <b>hadm[n+1][N]</b> 
116      * <br><i>(If I use [x][y] it always means [row-offset][column-offset] of a Matrix with n rows and m columns. Its entries go from M[0][0] to M[n][m])</i></li>
117      * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li>
118      * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows.
119      * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1].
120      * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column 
121      * </li>
122      * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows.
123      * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1].
124      * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column 
125      * </li>
126      * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above.
127      * <li>The output vector y is now in the last column of <b>hadm</b></li>
128      * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li>    
129      * </ol>
130      * <br>
131      * <b><u>Visually</u></b>
132      * <pre>
133      *        +--------+---+---+---+-----+---+
134      *        |   0    | 1 | 2 | 3 | ... |n+1|
135      * +------+--------+---+---+---+-----+---+
136      * |0     | x<sub>0</sub>     |       /\            |
137      * |1     | x<sub>1</sub>     |       ||            |
138      * |2     | x<sub>2</sub>     |   <= D<sub>top</sub>  =>       |
139      * |...   | ...    |       ||            |
140      * |N/2-1 | x<sub>N/2-1</sub>  |       \/            |
141      * +------+--------+---+---+---+-----+---+
142      * |N/2   | x<sub>N/2</sub>   |       /\            |
143      * |N/2+1 | x<sub>N/2+1</sub>  |       ||            |
144      * |N/2+2 | x<sub>N/2+2</sub>  |  <= D<sub>bottom</sub>  =>      | which is in the last column of the matrix
145      * |...   | ...    |       ||            |
146      * |N     | x<sub>N/2</sub>   |        \/           |
147      * +------+--------+---+---+---+-----+---+
148      * </pre>
149      * 
150      * @param x input vector
151      * @return y output vector
152      * @exception IllegalArgumentException if input array is not a power of 2
153      */
154     protected double[] fht(double x[]) throws IllegalArgumentException {
155 
156         // n is the row count of the input vector x
157         final int n     = x.length;
158         final int halfN = n / 2;
159 
160         // n has to be of the form n = 2^p !!
161         if (!FastFourierTransformer.isPowerOf2(n)) {
162             throw MathRuntimeException.createIllegalArgumentException(
163                     "{0} is not a power of 2",
164                     n);
165         }
166 
167         // Instead of creating a matrix with p+1 columns and n rows
168         // we will use two single dimension arrays which we will use in an alternating way.
169         double[] yPrevious = new double[n];
170         double[] yCurrent  = x.clone();
171 
172         // iterate from left to right (column)
173         for (int j = 1; j < n; j <<= 1) {
174 
175             // switch columns
176             final double[] yTmp = yCurrent;
177             yCurrent  = yPrevious;
178             yPrevious = yTmp;
179 
180             // iterate from top to bottom (row)
181             for (int i = 0; i < halfN; ++i) { 
182                 // D<sub>top</sub>
183                 // The top part works with addition
184                 final int twoI = 2 * i;
185                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
186             }
187             for (int i = halfN; i < n; ++i) { 
188                 // D<sub>bottom</sub>   
189                 // The bottom part works with subtraction
190                 final int twoI = 2 * i;
191                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
192             }
193         }
194 
195         // return the last computed output vector y
196         return yCurrent;
197 
198     }
199     /**
200      * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
201      * @param x input vector
202      * @return y output vector
203      * @exception IllegalArgumentException if input array is not a power of 2
204      */
205     protected int[] fht(int x[]) throws IllegalArgumentException {
206 
207         // n is the row count of the input vector x
208         final int n     = x.length;
209         final int halfN = n / 2;
210 
211         // n has to be of the form n = 2^p !!
212         if (!FastFourierTransformer.isPowerOf2(n)) {
213             throw MathRuntimeException.createIllegalArgumentException(
214                     "{0} is not a power of 2",
215                     n);
216         }
217 
218         // Instead of creating a matrix with p+1 columns and n rows
219         // we will use two single dimension arrays which we will use in an alternating way.
220         int[] yPrevious = new int[n];
221         int[] yCurrent  = x.clone();
222 
223         // iterate from left to right (column)
224         for (int j = 1; j < n; j <<= 1) {
225 
226             // switch columns
227             final int[] yTmp = yCurrent;
228             yCurrent  = yPrevious;
229             yPrevious = yTmp;
230 
231             // iterate from top to bottom (row)
232             for (int i = 0; i < halfN; ++i) { 
233                 // D<sub>top</sub>
234                 // The top part works with addition
235                 final int twoI = 2 * i;
236                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
237             }
238             for (int i = halfN; i < n; ++i) { 
239                 // D<sub>bottom</sub>   
240                 // The bottom part works with subtraction
241                 final int twoI = 2 * i;
242                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
243             }
244         }
245 
246         // return the last computed output vector y
247         return yCurrent;
248 
249     }
250 
251 }