001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.transform;
018    
019    import org.apache.commons.math.FunctionEvaluationException;
020    import org.apache.commons.math.MathRuntimeException;
021    import org.apache.commons.math.analysis.UnivariateRealFunction;
022    
023    /**
024     * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
025     * Transformation of an input vector x to the output vector y.
026     * <p>In addition to transformation of real vectors, the Hadamard transform can
027     * transform integer vectors into integer vectors. However, this integer transform
028     * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
029     * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
030     * vector (1/2, -1/2, 0, 0).</p>
031     * @version $Revision: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 2009) $
032     * @since 2.0
033     */
034    public class FastHadamardTransformer implements RealTransformer {
035    
036        /** {@inheritDoc} */
037        public double[] transform(double f[])
038            throws IllegalArgumentException {
039            return fht(f);
040        }
041    
042        /** {@inheritDoc} */
043        public double[] transform(UnivariateRealFunction f,
044                                  double min, double max, int n)
045            throws FunctionEvaluationException, IllegalArgumentException {
046            return fht(FastFourierTransformer.sample(f, min, max, n));
047        }
048    
049        /** {@inheritDoc} */
050        public double[] inversetransform(double f[])
051        throws IllegalArgumentException {
052            return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
053       }
054    
055        /** {@inheritDoc} */
056        public double[] inversetransform(UnivariateRealFunction f,
057                                         double min, double max, int n)
058            throws FunctionEvaluationException, IllegalArgumentException {
059            final double[] unscaled =
060                fht(FastFourierTransformer.sample(f, min, max, n));
061            return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
062        }
063    
064        /**
065         * Transform the given real data set.
066         * <p>The integer transform cannot be inverted directly, due to a scaling
067         * factor it may lead to double results.</p>
068         * @param f the integer data array to be transformed (signal)
069         * @return the integer transformed array (spectrum)
070         * @throws IllegalArgumentException if any parameters are invalid
071         */
072        public int[] transform(int f[])
073            throws IllegalArgumentException {
074            return fht(f);
075        }
076    
077        /**
078         * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
079         * <br>
080         * Requires <b>Nlog2N = n2</b><sup>n</sup> additions.
081         * <br>
082         * <br>
083         * <b><u>Short Table of manual calculation for N=8:</u></b>
084         * <ol>
085         * <li><b>x</b> is the input vector we want to transform</li>
086         * <li><b>y</b> is the output vector which is our desired result</li>
087         * <li>a and b are just helper rows</li>
088         * </ol>
089         * <pre>
090         * <code>
091         * +----+----------+---------+----------+
092         * | <b>x</b>  |    <b>a</b>     |    <b>b</b>    |    <b>y</b>     |
093         * +----+----------+---------+----------+
094         * | 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> |
095         * +----+----------+---------+----------+
096         * | 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> |
097         * +----+----------+---------+----------+
098         * | 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> |
099         * +----+----------+---------+----------+
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    }