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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.transform;

  24. import java.io.Serializable;

  26. import org.hipparchus.analysis.FunctionUtils;
  27. import org.hipparchus.analysis.UnivariateFunction;
  28. import org.hipparchus.exception.MathIllegalArgumentException;
  29. import org.hipparchus.util.ArithmeticUtils;

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

  44.     /** Serializable version identifier. */
  45.     static final long serialVersionUID = 20120211L;

  47.     /** Empty constructor.
  48.      * <p>
  49.      * This constructor is not strictly necessary, but it prevents spurious
  50.      * javadoc warnings with JDK 18 and later.
  51.      * </p>
  52.      * @since 3.0
  53.      */
  54.     public FastHadamardTransformer() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
  55.         // nothing to do
  56.     }

  58.     /**
  59.      * {@inheritDoc}
  60.      *
  61.      * @throws MathIllegalArgumentException if the length of the data array is
  62.      * not a power of two
  63.      */
  64.     @Override
  65.     public double[] transform(final double[] f, final TransformType type) {
  66.         if (type == TransformType.FORWARD) {
  67.             return fht(f);
  68.         }
  69.         return TransformUtils.scaleArray(fht(f), 1.0 / f.length);
  70.     }

  72.     /**
  73.      * {@inheritDoc}
  74.      *
  75.      * @throws org.hipparchus.exception.MathIllegalArgumentException
  76.      *   if the lower bound is greater than, or equal to the upper bound
  77.      * @throws org.hipparchus.exception.MathIllegalArgumentException
  78.      *   if the number of sample points is negative
  79.      * @throws MathIllegalArgumentException if the number of sample points is not a power of two
  80.      */
  81.     @Override
  82.     public double[] transform(final UnivariateFunction f,
  83.         final double min, final double max, final int n,
  84.         final TransformType type) {

  86.         return transform(FunctionUtils.sample(f, min, max, n), type);
  87.     }

  89.     /**
  90.      * Returns the forward transform of the specified integer data set.The
  91.      * integer transform cannot be inverted directly, due to a scaling factor
  92.      * which may lead to double results.
  93.      *
  94.      * @param f the integer data array to be transformed (signal)
  95.      * @return the integer transformed array (spectrum)
  96.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  97.      */
  98.     public int[] transform(final int[] f) {
  99.         return fht(f);
  100.     }

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

  246.         final int n     = x.length;
  247.         final int halfN = n / 2;

  249.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  250.             throw new MathIllegalArgumentException(LocalizedFFTFormats.NOT_POWER_OF_TWO,
  251.                     n);
  252.         }

  254.         /*
  255.          * Instead of creating a matrix with p+1 columns and n rows, we use two
  256.          * one dimension arrays which we are used in an alternating way.
  257.          */
  258.         double[] yPrevious = new double[n];
  259.         double[] yCurrent  = x.clone();

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

  264.             // switch columns
  265.             final double[] yTmp = yCurrent;
  266.             yCurrent  = yPrevious;
  267.             yPrevious = yTmp;

  269.             // iterate from top to bottom (row)
  270.             for (int i = 0; i < halfN; ++i) {
  271.                 // Dtop: the top part works with addition
  272.                 final int twoI = 2 * i;
  273.                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
  274.             }
  275.             for (int i = halfN; i < n; ++i) {
  276.                 // Dbottom: the bottom part works with subtraction
  277.                 final int twoI = 2 * i;
  278.                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
  279.             }
  280.         }

  282.         return yCurrent;

  284.     }

  286.     /**
  287.      * Returns the forward transform of the specified integer data set. The FHT
  288.      * (Fast Hadamard Transform) uses only subtraction and addition.
  289.      *
  290.      * @param x the integer data array to be transformed
  291.      * @return the integer transformed array, {@code y}
  292.      * @throws MathIllegalArgumentException if the length of the data array is not a power of two
  293.      */
  294.     protected int[] fht(int[] x) throws MathIllegalArgumentException {

  296.         final int n     = x.length;
  297.         final int halfN = n / 2;

  299.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  300.             throw new MathIllegalArgumentException(LocalizedFFTFormats.NOT_POWER_OF_TWO,
  301.                     n);
  302.         }

  304.         /*
  305.          * Instead of creating a matrix with p+1 columns and n rows, we use two
  306.          * one dimension arrays which we are used in an alternating way.
  307.          */
  308.         int[] yPrevious = new int[n];
  309.         int[] yCurrent  = x.clone();

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

  314.             // switch columns
  315.             final int[] yTmp = yCurrent;
  316.             yCurrent  = yPrevious;
  317.             yPrevious = yTmp;

  319.             // iterate from top to bottom (row)
  320.             for (int i = 0; i < halfN; ++i) {
  321.                 // Dtop: the top part works with addition
  322.                 final int twoI = 2 * i;
  323.                 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
  324.             }
  325.             for (int i = halfN; i < n; ++i) {
  326.                 // Dbottom: the bottom part works with subtraction
  327.                 final int twoI = 2 * i;
  328.                 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
  329.             }
  330.         }

  332.         // return the last computed output vector y
  333.         return yCurrent;

  335.     }

  337. }
Created with JaCoCo 0.8.12.202403310830