    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    import org.apache.commons.math.exception.util.LocalizedFormats;
    023
    024    /**
    025     * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
    026     * Transformation of an input vector x to the output vector y.
    027     * <p>In addition to transformation of real vectors, the Hadamard transform can
    028     * transform integer vectors into integer vectors. However, this integer transform
    029     * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
    030     * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
    031     * vector (1/2, -1/2, 0, 0).</p>
    032     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    033     * @since 2.0
    034     */
    035    public class FastHadamardTransformer implements RealTransformer {
    036
    037        /** {@inheritDoc} */
    038        public double[] transform(double f[])
    039            throws IllegalArgumentException {
    040            return fht(f);
    041        }
    042
    043        /** {@inheritDoc} */
    044        public double[] transform(UnivariateRealFunction f,
    045                                  double min, double max, int n)
    046            throws FunctionEvaluationException, IllegalArgumentException {
    047            return fht(FastFourierTransformer.sample(f, min, max, n));
    048        }
    049
    050        /** {@inheritDoc} */
    051        public double[] inversetransform(double f[])
    052        throws IllegalArgumentException {
    053            return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
    054       }
    055
    056        /** {@inheritDoc} */
    057        public double[] inversetransform(UnivariateRealFunction f,
    058                                         double min, double max, int n)
    059            throws FunctionEvaluationException, IllegalArgumentException {
    060            final double[] unscaled =
    061                fht(FastFourierTransformer.sample(f, min, max, n));
    062            return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
    063        }
    064
    065        /**
    066         * Transform the given real data set.
    067         * <p>The integer transform cannot be inverted directly, due to a scaling
    068         * factor it may lead to double results.</p>
    069         * @param f the integer data array to be transformed (signal)
    070         * @return the integer transformed array (spectrum)
    071         * @throws IllegalArgumentException if any parameters are invalid
    072         */
    073        public int[] transform(int f[])
    074            throws IllegalArgumentException {
    075            return fht(f);
    076        }
    077
    078        /**
    079         * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
    080         * <br>
    081         * Requires <b>Nlog2N = n2</b><sup>n</sup> additions.
    082         * <br>
    083         * <br>
    084         * <b><u>Short Table of manual calculation for N=8:</u></b>
    085         * <ol>
    086         * <li><b>x</b> is the input vector we want to transform</li>
    087         * <li><b>y</b> is the output vector which is our desired result</li>
    088         * <li>a and b are just helper rows</li>
    089         * </ol>
    090         * <pre>
    091         * <code>
    092         * +----+----------+---------+----------+
    093         * | <b>x</b>  |    <b>a</b>     |    <b>b</b>    |    <b>y</b>     |
    094         * +----+----------+---------+----------+
    095         * | 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> |
    096         * +----+----------+---------+----------+
    097         * | 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> |
    098         * +----+----------+---------+----------+
    099         * | 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> |
    100         * +----+----------+---------+----------+
    101         * | 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> |
    102         * +----+----------+---------+----------+
    103         * | 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> |
    104         * +----+----------+---------+----------+
    105         * | 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> |
    106         * +----+----------+---------+----------+
    107         * | 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> |
    108         * +----+----------+---------+----------+
    109         * | 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> |
    110         * +----+----------+---------+----------+
    111         * </code>
    112         * </pre>
    113         *
    114         * <b><u>How it works</u></b>
    115         * <ol>
    116         * <li>Construct a matrix with N rows and n+1 columns<br>   <b>hadm[n+1][N]</b>
    117         * <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>
    118         * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li>
    119         * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows.
    120         * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1].
    121         * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column
    122         * </li>
    123         * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows.
    124         * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1].
    125         * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column
    126         * </li>
    127         * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above.
    128         * <li>The output vector y is now in the last column of <b>hadm</b></li>
    129         * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li>
    130         * </ol>
    131         * <br>
    132         * <b><u>Visually</u></b>
    133         * <pre>
    134         *        +--------+---+---+---+-----+---+
    135         *        |   0    | 1 | 2 | 3 | ... |n+1|
    136         * +------+--------+---+---+---+-----+---+
    137         * |0     | x<sub>0</sub>     |       /\            |
    138         * |1     | x<sub>1</sub>     |       ||            |
    139         * |2     | x<sub>2</sub>     |   <= D<sub>top</sub>  =>       |
    140         * |...   | ...    |       ||            |
    141         * |N/2-1 | x<sub>N/2-1</sub>  |       \/            |
    142         * +------+--------+---+---+---+-----+---+
    143         * |N/2   | x<sub>N/2</sub>   |       /\            |
    144         * |N/2+1 | x<sub>N/2+1</sub>  |       ||            |
    145         * |N/2+2 | x<sub>N/2+2</sub>  |  <= D<sub>bottom</sub>  =>      | which is in the last column of the matrix
    146         * |...   | ...    |       ||            |
    147         * |N     | x<sub>N/2</sub>   |        \/           |
    148         * +------+--------+---+---+---+-----+---+
    149         * </pre>
    150         *
    151         * @param x input vector
    152         * @return y output vector
    153         * @exception IllegalArgumentException if input array is not a power of 2
    154         */
    155        protected double[] fht(double x[]) throws IllegalArgumentException {
    156
    157            // n is the row count of the input vector x
    158            final int n     = x.length;
    159            final int halfN = n / 2;
    160
    161            // n has to be of the form n = 2^p !!
    162            if (!FastFourierTransformer.isPowerOf2(n)) {
    163                throw MathRuntimeException.createIllegalArgumentException(
    164                        LocalizedFormats.NOT_POWER_OF_TWO,
    165                        n);
    166            }
    167
    168            // Instead of creating a matrix with p+1 columns and n rows
    169            // we will use two single dimension arrays which we will use in an alternating way.
    170            double[] yPrevious = new double[n];
    171            double[] yCurrent  = x.clone();
    172
    173            // iterate from left to right (column)
    174            for (int j = 1; j < n; j <<= 1) {
    175
    176                // switch columns
    177                final double[] yTmp = yCurrent;
    178                yCurrent  = yPrevious;
    179                yPrevious = yTmp;
    180
    181                // iterate from top to bottom (row)
    182                for (int i = 0; i < halfN; ++i) {
    183                    // D<sub>top</sub>
    184                    // The top part works with addition
    185                    final int twoI = 2 * i;
    186                    yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
    187                }
    188                for (int i = halfN; i < n; ++i) {
    189                    // D<sub>bottom</sub>
    190                    // The bottom part works with subtraction
    191                    final int twoI = 2 * i;
    192                    yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
    193                }
    194            }
    195
    196            // return the last computed output vector y
    197            return yCurrent;
    198
    199        }
    200        /**
    201         * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
    202         * @param x input vector
    203         * @return y output vector
    204         * @exception IllegalArgumentException if input array is not a power of 2
    205         */
    206        protected int[] fht(int x[]) throws IllegalArgumentException {
    207
    208            // n is the row count of the input vector x
    209            final int n     = x.length;
    210            final int halfN = n / 2;
    211
    212            // n has to be of the form n = 2^p !!
    213            if (!FastFourierTransformer.isPowerOf2(n)) {
    214                throw MathRuntimeException.createIllegalArgumentException(
    215                        LocalizedFormats.NOT_POWER_OF_TWO,
    216                        n);
    217            }
    218
    219            // Instead of creating a matrix with p+1 columns and n rows
    220            // we will use two single dimension arrays which we will use in an alternating way.
    221            int[] yPrevious = new int[n];
    222            int[] yCurrent  = x.clone();
    223
    224            // iterate from left to right (column)
    225            for (int j = 1; j < n; j <<= 1) {
    226
    227                // switch columns
    228                final int[] yTmp = yCurrent;
    229                yCurrent  = yPrevious;
    230                yPrevious = yTmp;
    231
    232                // iterate from top to bottom (row)
    233                for (int i = 0; i < halfN; ++i) {
    234                    // D<sub>top</sub>
    235                    // The top part works with addition
    236                    final int twoI = 2 * i;
    237                    yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
    238                }
    239                for (int i = halfN; i < n; ++i) {
    240                    // D<sub>bottom</sub>
    241                    // The bottom part works with subtraction
    242                    final int twoI = 2 * i;
    243                    yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
    244                }
    245            }
    246
    247            // return the last computed output vector y
    248            return yCurrent;
    249
    250        }
    251
    252    }