FastHadamardTransformer.java

  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.math3.transform;

  18. import java.io.Serializable;

  19. import org.apache.commons.math3.analysis.FunctionUtils;
  20. import org.apache.commons.math3.analysis.UnivariateFunction;
  21. import org.apache.commons.math3.exception.MathIllegalArgumentException;
  22. import org.apache.commons.math3.exception.util.LocalizedFormats;
  23. import org.apache.commons.math3.util.ArithmeticUtils;

  24. /**
  25.  * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
  26.  * Transformation of an input vector x to the output vector y.
  27.  * <p>
  28.  * In addition to transformation of real vectors, the Hadamard transform can
  29.  * transform integer vectors into integer vectors. However, this integer transform
  30.  * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
  31.  * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
  32.  * vector (1/2, -1/2, 0, 0).
  33.  *
  34.  * @since 2.0
  35.  */
  36. public class FastHadamardTransformer implements RealTransformer, Serializable {

  37.     /** Serializable version identifier. */
  38.     static final long serialVersionUID = 20120211L;

  39.     /**
  40.      * {@inheritDoc}
  41.      *
  42.      * @throws MathIllegalArgumentException if the length of the data array is
  43.      * not a power of two
  44.      */
  45.     public double[] transform(final double[] f, final TransformType type) {
  46.         if (type == TransformType.FORWARD) {
  47.             return fht(f);
  48.         }
  49.         return TransformUtils.scaleArray(fht(f), 1.0 / f.length);
  50.     }

  51.     /**
  52.      * {@inheritDoc}
  53.      *
  54.      * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
  55.      *   if the lower bound is greater than, or equal to the upper bound
  56.      * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
  57.      *   if the number of sample points is negative
  58.      * @throws MathIllegalArgumentException if the number of sample points is not a power of two
  59.      */
  60.     public double[] transform(final UnivariateFunction f,
  61.         final double min, final double max, final int n,
  62.         final TransformType type) {

  63.         return transform(FunctionUtils.sample(f, min, max, n), type);
  64.     }

  65.     /**
  66.      * Returns the forward transform of the specified integer data set.The
  67.      * integer transform cannot be inverted directly, due to a scaling factor
  68.      * which may lead to double results.
  69.      *
  70.      * @param f the integer data array to be transformed (signal)
  71.      * @return the integer transformed array (spectrum)
  72.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  73.      */
  74.     public int[] transform(final int[] f) {
  75.         return fht(f);
  76.     }

  77.     /**
  78.      * The FHT (Fast Hadamard Transformation) which uses only subtraction and
  79.      * addition. Requires {@code N * log2(N) = n * 2^n} additions.
  80.      *
  81.      * <h3>Short Table of manual calculation for N=8</h3>
  82.      * <ol>
  83.      * <li><b>x</b> is the input vector to be transformed,</li>
  84.      * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>),</li>
  85.      * <li>a and b are helper rows.</li>
  86.      * </ol>
  87.      * <table align="center" border="1" cellpadding="3">
  88.      * <tbody align="center">
  89.      * <tr>
  90.      *     <th>x</th>
  91.      *     <th>a</th>
  92.      *     <th>b</th>
  93.      *     <th>y</th>
  94.      * </tr>
  95.      * <tr>
  96.      *     <th>x<sub>0</sub></th>
  97.      *     <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td>
  98.      *     <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td>
  99.      *     <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td>
  100.      * </tr>
  101.      * <tr>
  102.      *     <th>x<sub>1</sub></th>
  103.      *     <td>a<sub>1</sub> = x<sub>2</sub> + x<sub>3</sub></td>
  104.      *     <td>b<sub>0</sub> = a<sub>2</sub> + a<sub>3</sub></td>
  105.      *     <td>y<sub>0</sub> = b<sub>2</sub> + b<sub>3</sub></td>
  106.      * </tr>
  107.      * <tr>
  108.      *     <th>x<sub>2</sub></th>
  109.      *     <td>a<sub>2</sub> = x<sub>4</sub> + x<sub>5</sub></td>
  110.      *     <td>b<sub>0</sub> = a<sub>4</sub> + a<sub>5</sub></td>
  111.      *     <td>y<sub>0</sub> = b<sub>4</sub> + b<sub>5</sub></td>
  112.      * </tr>
  113.      * <tr>
  114.      *     <th>x<sub>3</sub></th>
  115.      *     <td>a<sub>3</sub> = x<sub>6</sub> + x<sub>7</sub></td>
  116.      *     <td>b<sub>0</sub> = a<sub>6</sub> + a<sub>7</sub></td>
  117.      *     <td>y<sub>0</sub> = b<sub>6</sub> + b<sub>7</sub></td>
  118.      * </tr>
  119.      * <tr>
  120.      *     <th>x<sub>4</sub></th>
  121.      *     <td>a<sub>0</sub> = x<sub>0</sub> - x<sub>1</sub></td>
  122.      *     <td>b<sub>0</sub> = a<sub>0</sub> - a<sub>1</sub></td>
  123.      *     <td>y<sub>0</sub> = b<sub>0</sub> - b<sub>1</sub></td>
  124.      * </tr>
  125.      * <tr>
  126.      *     <th>x<sub>5</sub></th>
  127.      *     <td>a<sub>1</sub> = x<sub>2</sub> - x<sub>3</sub></td>
  128.      *     <td>b<sub>0</sub> = a<sub>2</sub> - a<sub>3</sub></td>
  129.      *     <td>y<sub>0</sub> = b<sub>2</sub> - b<sub>3</sub></td>
  130.      * </tr>
  131.      * <tr>
  132.      *     <th>x<sub>6</sub></th>
  133.      *     <td>a<sub>2</sub> = x<sub>4</sub> - x<sub>5</sub></td>
  134.      *     <td>b<sub>0</sub> = a<sub>4</sub> - a<sub>5</sub></td>
  135.      *     <td>y<sub>0</sub> = b<sub>4</sub> - b<sub>5</sub></td>
  136.      * </tr>
  137.      * <tr>
  138.      *     <th>x<sub>7</sub></th>
  139.      *     <td>a<sub>3</sub> = x<sub>6</sub> - x<sub>7</sub></td>
  140.      *     <td>b<sub>0</sub> = a<sub>6</sub> - a<sub>7</sub></td>
  141.      *     <td>y<sub>0</sub> = b<sub>6</sub> - b<sub>7</sub></td>
  142.      * </tr>
  143.      * </tbody>
  144.      * </table>
  145.      *
  146.      * <h3>How it works</h3>
  147.      * <ol>
  148.      * <li>Construct a matrix with {@code N} rows and {@code n + 1} columns,
  149.      * {@code hadm[n+1][N]}.<br/>
  150.      * <em>(If I use [x][y] it always means [row-offset][column-offset] of a
  151.      * Matrix with n rows and m columns. Its entries go from M[0][0]
  152.      * to M[n][N])</em></li>
  153.      * <li>Place the input vector {@code x[N]} in the first column of the
  154.      * matrix {@code hadm}.</li>
  155.      * <li>The entries of the submatrix {@code D_top} are calculated as follows
  156.      *     <ul>
  157.      *         <li>{@code D_top} goes from entry {@code [0][1]} to
  158.      *         {@code [N / 2 - 1][n + 1]},</li>
  159.      *         <li>the columns of {@code D_top} are the pairwise mutually
  160.      *         exclusive sums of the previous column.</li>
  161.      *     </ul>
  162.      * </li>
  163.      * <li>The entries of the submatrix {@code D_bottom} are calculated as
  164.      * follows
  165.      *     <ul>
  166.      *         <li>{@code D_bottom} goes from entry {@code [N / 2][1]} to
  167.      *         {@code [N][n + 1]},</li>
  168.      *         <li>the columns of {@code D_bottom} are the pairwise differences
  169.      *         of the previous column.</li>
  170.      *     </ul>
  171.      * </li>
  172.      * <li>The consputation of {@code D_top} and {@code D_bottom} are best
  173.      * understood with the above example (for {@code N = 8}).
  174.      * <li>The output vector {@code y} is now in the last column of
  175.      * {@code hadm}.</li>
  176.      * <li><em>Algorithm from <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">chipcenter</a>.</em></li>
  177.      * </ol>
  178.      * <h3>Visually</h3>
  179.      * <table border="1" align="center" cellpadding="3">
  180.      * <tbody align="center">
  181.      * <tr>
  182.      *     <td></td><th>0</th><th>1</th><th>2</th><th>3</th>
  183.      *     <th>&hellip;</th>
  184.      *     <th>n + 1</th>
  185.      * </tr>
  186.      * <tr>
  187.      *     <th>0</th>
  188.      *     <td>x<sub>0</sub></td>
  189.      *     <td colspan="5" rowspan="5" align="center" valign="middle">
  190.      *         &uarr;<br/>
  191.      *         &larr; D<sub>top</sub> &rarr;<br/>
  192.      *         &darr;
  193.      *     </td>
  194.      * </tr>
  195.      * <tr><th>1</th><td>x<sub>1</sub></td></tr>
  196.      * <tr><th>2</th><td>x<sub>2</sub></td></tr>
  197.      * <tr><th>&hellip;</th><td>&hellip;</td></tr>
  198.      * <tr><th>N / 2 - 1</th><td>x<sub>N/2-1</sub></td></tr>
  199.      * <tr>
  200.      *     <th>N / 2</th>
  201.      *     <td>x<sub>N/2</sub></td>
  202.      *     <td colspan="5" rowspan="5" align="center" valign="middle">
  203.      *         &uarr;<br/>
  204.      *         &larr; D<sub>bottom</sub> &rarr;<br/>
  205.      *         &darr;
  206.      *     </td>
  207.      * </tr>
  208.      * <tr><th>N / 2 + 1</th><td>x<sub>N/2+1</sub></td></tr>
  209.      * <tr><th>N / 2 + 2</th><td>x<sub>N/2+2</sub></td></tr>
  210.      * <tr><th>&hellip;</th><td>&hellip;</td></tr>
  211.      * <tr><th>N</th><td>x<sub>N</sub></td></tr>
  212.      * </tbody>
  213.      * </table>
  214.      *
  215.      * @param x the real data array to be transformed
  216.      * @return the real transformed array, {@code y}
  217.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  218.      */
  219.     protected double[] fht(double[] x) throws MathIllegalArgumentException {

  220.         final int n     = x.length;
  221.         final int halfN = n / 2;

  222.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  223.             throw new MathIllegalArgumentException(
  224.                     LocalizedFormats.NOT_POWER_OF_TWO,
  225.                     Integer.valueOf(n));
  226.         }

  227.         /*
  228.          * Instead of creating a matrix with p+1 columns and n rows, we use two
  229.          * one dimension arrays which we are used in an alternating way.
  230.          */
  231.         double[] yPrevious = new double[n];
  232.         double[] yCurrent  = x.clone();

  233.         // iterate from left to right (column)
  234.         for (int j = 1; j < n; j <<= 1) {

  235.             // switch columns
  236.             final double[] yTmp = yCurrent;
  237.             yCurrent  = yPrevious;
  238.             yPrevious = yTmp;

  239.             // iterate from top to bottom (row)
  240.             for (int i = 0; i < halfN; ++i) {
  241.                 // Dtop: the top part works with addition
  242.                 final int twoI = 2 * i;
  243.                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
  244.             }
  245.             for (int i = halfN; i < n; ++i) {
  246.                 // Dbottom: the bottom part works with subtraction
  247.                 final int twoI = 2 * i;
  248.                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
  249.             }
  250.         }

  251.         return yCurrent;

  252.     }

  253.     /**
  254.      * Returns the forward transform of the specified integer data set. The FHT
  255.      * (Fast Hadamard Transform) uses only subtraction and addition.
  256.      *
  257.      * @param x the integer data array to be transformed
  258.      * @return the integer transformed array, {@code y}
  259.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  260.      */
  261.     protected int[] fht(int[] x) throws MathIllegalArgumentException {

  262.         final int n     = x.length;
  263.         final int halfN = n / 2;

  264.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  265.             throw new MathIllegalArgumentException(
  266.                     LocalizedFormats.NOT_POWER_OF_TWO,
  267.                     Integer.valueOf(n));
  268.         }

  269.         /*
  270.          * Instead of creating a matrix with p+1 columns and n rows, we use two
  271.          * one dimension arrays which we are used in an alternating way.
  272.          */
  273.         int[] yPrevious = new int[n];
  274.         int[] yCurrent  = x.clone();

  275.         // iterate from left to right (column)
  276.         for (int j = 1; j < n; j <<= 1) {

  277.             // switch columns
  278.             final int[] yTmp = yCurrent;
  279.             yCurrent  = yPrevious;
  280.             yPrevious = yTmp;

  281.             // iterate from top to bottom (row)
  282.             for (int i = 0; i < halfN; ++i) {
  283.                 // Dtop: the top part works with addition
  284.                 final int twoI = 2 * i;
  285.                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
  286.             }
  287.             for (int i = halfN; i < n; ++i) {
  288.                 // Dbottom: the bottom part works with subtraction
  289.                 final int twoI = 2 * i;
  290.                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
  291.             }
  292.         }

  293.         // return the last computed output vector y
  294.         return yCurrent;

  295.     }

  296. }