FastHadamardTransform.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.math4.transform;

  19. import java.util.function.UnaryOperator;

  21. import org.apache.commons.numbers.core.ArithmeticUtils;

  23. /**
  24.  * <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
  25.  * <p>
  26.  * The FHT can also transform integer vectors into integer vectors.
  27.  * However, this transform cannot be inverted directly, due to a scaling
  28.  * factor that may lead to rational results (for example, the inverse
  29.  * transform of integer vector (0, 1, 0, 1) is vector (1/2, -1/2, 0, 0).
  30.  */
  31. public class FastHadamardTransform implements RealTransform {
  32.     /** Operation to be performed. */
  33.     private final UnaryOperator<double[]> op;

  35.     /**
  36.      * Default constructor.
  37.      */
  38.     public FastHadamardTransform() {
  39.         this(false);
  40.     }

  42.     /**
  43.      * @param inverse Whether to perform the inverse transform.
  44.      */
  45.     public FastHadamardTransform(final boolean inverse) {
  46.         op = create(inverse);
  47.     }

  49.     /**
  50.      * {@inheritDoc}
  51.      *
  52.      * @throws IllegalArgumentException if the length of the data array is
  53.      * not a power of two.
  54.      */
  55.     @Override
  56.     public double[] apply(final double[] f) {
  57.         return op.apply(f);
  58.     }

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

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

  222.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  223.             throw new TransformException(TransformException.NOT_POWER_OF_TWO,
  224.                                          n);
  225.         }

  227.         final int halfN = n / 2;

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

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

  237.             // switch columns
  238.             final double[] yTmp = yCurrent;
  239.             yCurrent  = yPrevious;
  240.             yPrevious = yTmp;

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

  255.         return yCurrent;
  256.     }

  258.     /**
  259.      * Returns the forward transform of the specified integer data set.
  260.      * FHT only uses subtraction and addition.
  261.      *
  262.      * @param x Data to be transformed.
  263.      * @return the transformed array.
  264.      * @throws IllegalArgumentException if the length of the data array is
  265.      * not a power of two.
  266.      */
  267.     private int[] fht(int[] x) {
  268.         final int n = x.length;
  269.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  270.             throw new TransformException(TransformException.NOT_POWER_OF_TWO,
  271.                                          n);
  272.         }

  274.         final int halfN = n / 2;

  276.         // Instead of creating a matrix with p+1 columns and n rows, we use two
  277.         // one dimension arrays which we are used in an alternating way.

  279.         int[] yPrevious = new int[n];
  280.         int[] yCurrent  = x.clone();

  282.         // iterate from left to right (column)
  283.         for (int j = 1; j < n; j <<= 1) {
  284.             // switch columns
  285.             final int[] yTmp = yCurrent;
  286.             yCurrent  = yPrevious;
  287.             yPrevious = yTmp;

  289.             // iterate from top to bottom (row)
  290.             for (int i = 0; i < halfN; ++i) {
  291.                 // Dtop: the top part works with addition
  292.                 final int twoI = 2 * i;
  293.                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
  294.             }
  295.             for (int i = halfN; i < n; ++i) {
  296.                 // Dbottom: the bottom part works with subtraction
  297.                 final int twoI = 2 * i;
  298.                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
  299.             }
  300.         }

  302.         // return the last computed output vector y
  303.         return yCurrent;
  304.     }

  306.     /**
  307.      * Factory method.
  308.      *
  309.      * @param inverse Whether to perform the inverse transform.
  310.      * @return the transform operator.
  311.      */
  312.     private UnaryOperator<double[]> create(final boolean inverse) {
  313.         if (inverse) {
  314.             return f -> TransformUtils.scaleInPlace(fht(f), 1d / f.length);
  315.         } else {
  316.             return f -> fht(f);
  317.         }
  318.     }
  319. }
Created with JaCoCo 0.8.10.202304240956