    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.math3.transform;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math3.analysis.FunctionUtils;
    022    import org.apache.commons.math3.analysis.UnivariateFunction;
    023    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025    import org.apache.commons.math3.util.ArithmeticUtils;
    026
    027    /**
    028     * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
    029     * Transformation of an input vector x to the output vector y.
    030     * <p>
    031     * In addition to transformation of real vectors, the Hadamard transform can
    032     * transform integer vectors into integer vectors. However, this integer transform
    033     * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
    034     * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
    035     * vector (1/2, -1/2, 0, 0).
    036     *
    037     * @version $Id: FastHadamardTransformer.java 1385310 2012-09-16 16:32:10Z tn $
    038     * @since 2.0
    039     */
    040    public class FastHadamardTransformer implements RealTransformer, Serializable {
    041
    042        /** Serializable version identifier. */
    043        static final long serialVersionUID = 20120211L;
    044
    045        /**
    046         * {@inheritDoc}
    047         *
    048         * @throws MathIllegalArgumentException if the length of the data array is
    049         * not a power of two
    050         */
    051        public double[] transform(final double[] f, final TransformType type) {
    052            if (type == TransformType.FORWARD) {
    053                return fht(f);
    054            }
    055            return TransformUtils.scaleArray(fht(f), 1.0 / f.length);
    056        }
    057
    058        /**
    059         * {@inheritDoc}
    060         *
    061         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    062         *   if the lower bound is greater than, or equal to the upper bound
    063         * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
    064         *   if the number of sample points is negative
    065         * @throws MathIllegalArgumentException if the number of sample points is not a power of two
    066         */
    067        public double[] transform(final UnivariateFunction f,
    068            final double min, final double max, final int n,
    069            final TransformType type) {
    070
    071            return transform(FunctionUtils.sample(f, min, max, n), type);
    072        }
    073
    074        /**
    075         * Returns the forward transform of the specified integer data set.The
    076         * integer transform cannot be inverted directly, due to a scaling factor
    077         * which may lead to double results.
    078         *
    079         * @param f the integer data array to be transformed (signal)
    080         * @return the integer transformed array (spectrum)
    081         * @throws MathIllegalArgumentException if the length of the data array is not a power of two
    082         */
    083        public int[] transform(final int[] f) {
    084            return fht(f);
    085        }
    086
    087        /**
    088         * The FHT (Fast Hadamard Transformation) which uses only subtraction and
    089         * addition. Requires {@code N * log2(N) = n * 2^n} additions.
    090         *
    091         * <h3>Short Table of manual calculation for N=8</h3>
    092         * <ol>
    093         * <li><b>x</b> is the input vector to be transformed,</li>
    094         * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>),</li>
    095         * <li>a and b are helper rows.</li>
    096         * </ol>
    097         * <table align="center" border="1" cellpadding="3">
    098         * <tbody align="center">
    099         * <tr>
    100         *     <th>x</th>
    101         *     <th>a</th>
    102         *     <th>b</th>
    103         *     <th>y</th>
    104         * </tr>
    105         * <tr>
    106         *     <th>x<sub>0</sub></th>
    107         *     <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td>
    108         *     <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td>
    109         *     <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td>
    110         * </tr>
    111         * <tr>
    112         *     <th>x<sub>1</sub></th>
    113         *     <td>a<sub>1</sub> = x<sub>2</sub> + x<sub>3</sub></td>
    114         *     <td>b<sub>0</sub> = a<sub>2</sub> + a<sub>3</sub></td>
    115         *     <td>y<sub>0</sub> = b<sub>2</sub> + b<sub>3</sub></td>
    116         * </tr>
    117         * <tr>
    118         *     <th>x<sub>2</sub></th>
    119         *     <td>a<sub>2</sub> = x<sub>4</sub> + x<sub>5</sub></td>
    120         *     <td>b<sub>0</sub> = a<sub>4</sub> + a<sub>5</sub></td>
    121         *     <td>y<sub>0</sub> = b<sub>4</sub> + b<sub>5</sub></td>
    122         * </tr>
    123         * <tr>
    124         *     <th>x<sub>3</sub></th>
    125         *     <td>a<sub>3</sub> = x<sub>6</sub> + x<sub>7</sub></td>
    126         *     <td>b<sub>0</sub> = a<sub>6</sub> + a<sub>7</sub></td>
    127         *     <td>y<sub>0</sub> = b<sub>6</sub> + b<sub>7</sub></td>
    128         * </tr>
    129         * <tr>
    130         *     <th>x<sub>4</sub></th>
    131         *     <td>a<sub>0</sub> = x<sub>0</sub> - x<sub>1</sub></td>
    132         *     <td>b<sub>0</sub> = a<sub>0</sub> - a<sub>1</sub></td>
    133         *     <td>y<sub>0</sub> = b<sub>0</sub> - b<sub>1</sub></td>
    134         * </tr>
    135         * <tr>
    136         *     <th>x<sub>5</sub></th>
    137         *     <td>a<sub>1</sub> = x<sub>2</sub> - x<sub>3</sub></td>
    138         *     <td>b<sub>0</sub> = a<sub>2</sub> - a<sub>3</sub></td>
    139         *     <td>y<sub>0</sub> = b<sub>2</sub> - b<sub>3</sub></td>
    140         * </tr>
    141         * <tr>
    142         *     <th>x<sub>6</sub></th>
    143         *     <td>a<sub>2</sub> = x<sub>4</sub> - x<sub>5</sub></td>
    144         *     <td>b<sub>0</sub> = a<sub>4</sub> - a<sub>5</sub></td>
    145         *     <td>y<sub>0</sub> = b<sub>4</sub> - b<sub>5</sub></td>
    146         * </tr>
    147         * <tr>
    148         *     <th>x<sub>7</sub></th>
    149         *     <td>a<sub>3</sub> = x<sub>6</sub> - x<sub>7</sub></td>
    150         *     <td>b<sub>0</sub> = a<sub>6</sub> - a<sub>7</sub></td>
    151         *     <td>y<sub>0</sub> = b<sub>6</sub> - b<sub>7</sub></td>
    152         * </tr>
    153         * </tbody>
    154         * </table>
    155         *
    156         * <h3>How it works</h3>
    157         * <ol>
    158         * <li>Construct a matrix with {@code N} rows and {@code n + 1} columns,
    159         * {@code hadm[n+1][N]}.<br/>
    160         * <em>(If I use [x][y] it always means [row-offset][column-offset] of a
    161         * Matrix with n rows and m columns. Its entries go from M[0][0]
    162         * to M[n][N])</em></li>
    163         * <li>Place the input vector {@code x[N]} in the first column of the
    164         * matrix {@code hadm}.</li>
    165         * <li>The entries of the submatrix {@code D_top} are calculated as follows
    166         *     <ul>
    167         *         <li>{@code D_top} goes from entry {@code [0][1]} to
    168         *         {@code [N / 2 - 1][n + 1]},</li>
    169         *         <li>the columns of {@code D_top} are the pairwise mutually
    170         *         exclusive sums of the previous column.</li>
    171         *     </ul>
    172         * </li>
    173         * <li>The entries of the submatrix {@code D_bottom} are calculated as
    174         * follows
    175         *     <ul>
    176         *         <li>{@code D_bottom} goes from entry {@code [N / 2][1]} to
    177         *         {@code [N][n + 1]},</li>
    178         *         <li>the columns of {@code D_bottom} are the pairwise differences
    179         *         of the previous column.</li>
    180         *     </ul>
    181         * </li>
    182         * <li>The consputation of {@code D_top} and {@code D_bottom} are best
    183         * understood with the above example (for {@code N = 8}).
    184         * <li>The output vector {@code y} is now in the last column of
    185         * {@code hadm}.</li>
    186         * <li><em>Algorithm from <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">chipcenter</a>.</em></li>
    187         * </ol>
    188         * <h3>Visually</h3>
    189         * <table border="1" align="center" cellpadding="3">
    190         * <tbody align="center">
    191         * <tr>
    192         *     <td></td><th>0</th><th>1</th><th>2</th><th>3</th>
    193         *     <th>&hellip;</th>
    194         *     <th>n + 1</th>
    195         * </tr>
    196         * <tr>
    197         *     <th>0</th>
    198         *     <td>x<sub>0</sub></td>
    199         *     <td colspan="5" rowspan="5" align="center" valign="middle">
    200         *         &uarr;<br/>
    201         *         &larr; D<sub>top</sub> &rarr;<br/>
    202         *         &darr;
    203         *     </td>
    204         * </tr>
    205         * <tr><th>1</th><td>x<sub>1</sub></td></tr>
    206         * <tr><th>2</th><td>x<sub>2</sub></td></tr>
    207         * <tr><th>&hellip;</th><td>&hellip;</td></tr>
    208         * <tr><th>N / 2 - 1</th><td>x<sub>N/2-1</sub></td></tr>
    209         * <tr>
    210         *     <th>N / 2</th>
    211         *     <td>x<sub>N/2</sub></td>
    212         *     <td colspan="5" rowspan="5" align="center" valign="middle">
    213         *         &uarr;<br/>
    214         *         &larr; D<sub>bottom</sub> &rarr;<br/>
    215         *         &darr;
    216         *     </td>
    217         * </tr>
    218         * <tr><th>N / 2 + 1</th><td>x<sub>N/2+1</sub></td></tr>
    219         * <tr><th>N / 2 + 2</th><td>x<sub>N/2+2</sub></td></tr>
    220         * <tr><th>&hellip;</th><td>&hellip;</td></tr>
    221         * <tr><th>N</th><td>x<sub>N</sub></td></tr>
    222         * </tbody>
    223         * </table>
    224         *
    225         * @param x the real data array to be transformed
    226         * @return the real transformed array, {@code y}
    227         * @throws MathIllegalArgumentException if the length of the data array is not a power of two
    228         */
    229        protected double[] fht(double[] x) throws MathIllegalArgumentException {
    230
    231            final int n     = x.length;
    232            final int halfN = n / 2;
    233
    234            if (!ArithmeticUtils.isPowerOfTwo(n)) {
    235                throw new MathIllegalArgumentException(
    236                        LocalizedFormats.NOT_POWER_OF_TWO,
    237                        Integer.valueOf(n));
    238            }
    239
    240            /*
    241             * Instead of creating a matrix with p+1 columns and n rows, we use two
    242             * one dimension arrays which we are used in an alternating way.
    243             */
    244            double[] yPrevious = new double[n];
    245            double[] yCurrent  = x.clone();
    246
    247            // iterate from left to right (column)
    248            for (int j = 1; j < n; j <<= 1) {
    249
    250                // switch columns
    251                final double[] yTmp = yCurrent;
    252                yCurrent  = yPrevious;
    253                yPrevious = yTmp;
    254
    255                // iterate from top to bottom (row)
    256                for (int i = 0; i < halfN; ++i) {
    257                    // Dtop: the top part works with addition
    258                    final int twoI = 2 * i;
    259                    yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
    260                }
    261                for (int i = halfN; i < n; ++i) {
    262                    // Dbottom: the bottom part works with subtraction
    263                    final int twoI = 2 * i;
    264                    yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
    265                }
    266            }
    267
    268            return yCurrent;
    269
    270        }
    271
    272        /**
    273         * Returns the forward transform of the specified integer data set. The FHT
    274         * (Fast Hadamard Transform) uses only subtraction and addition.
    275         *
    276         * @param x the integer data array to be transformed
    277         * @return the integer transformed array, {@code y}
    278         * @throws MathIllegalArgumentException if the length of the data array is not a power of two
    279         */
    280        protected int[] fht(int[] x) throws MathIllegalArgumentException {
    281
    282            final int n     = x.length;
    283            final int halfN = n / 2;
    284
    285            if (!ArithmeticUtils.isPowerOfTwo(n)) {
    286                throw new MathIllegalArgumentException(
    287                        LocalizedFormats.NOT_POWER_OF_TWO,
    288                        Integer.valueOf(n));
    289            }
    290
    291            /*
    292             * Instead of creating a matrix with p+1 columns and n rows, we use two
    293             * one dimension arrays which we are used in an alternating way.
    294             */
    295            int[] yPrevious = new int[n];
    296            int[] yCurrent  = x.clone();
    297
    298            // iterate from left to right (column)
    299            for (int j = 1; j < n; j <<= 1) {
    300
    301                // switch columns
    302                final int[] yTmp = yCurrent;
    303                yCurrent  = yPrevious;
    304                yPrevious = yTmp;
    305
    306                // iterate from top to bottom (row)
    307                for (int i = 0; i < halfN; ++i) {
    308                    // Dtop: the top part works with addition
    309                    final int twoI = 2 * i;
    310                    yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
    311                }
    312                for (int i = halfN; i < n; ++i) {
    313                    // Dbottom: the bottom part works with subtraction
    314                    final int twoI = 2 * i;
    315                    yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
    316                }
    317            }
    318
    319            // return the last computed output vector y
    320            return yCurrent;
    321
    322        }
    323
    324    }