http://commons.apache.org/proper/commons-math/: The Apache Commons Math project is a library of lightweight, self-contained mathematics and statistics components addressing the most common practical problems not immediately available in the Java programming language or commons-lang. (The Apache Software Foundation)
  • Eldar Agalarov
  • Tim Allison
  • C. Scott Ananian
  • Mark Anderson
  • Peter Andrews
  • Rémi Arntzen
  • Matt Adereth
  • Jared Becksfort
  • Michael Bjorkegren
  • Brian Bloniarz
  • John Bollinger
  • Cyril Briquet
  • Dave Brosius
  • Dan Checkoway
  • Anders Conbere
  • Charles Cooper
  • Paul Cowan
  • Benjamin Croizet
  • Larry Diamond
  • Aleksei Dievskii
  • Rodrigo di Lorenzo Lopes
  • Hasan Diwan
  • Ted Dunning
  • Ole Ersoy
  • Ajo Fod
  • John Gant
  • Ken Geis
  • Hank Grabowski
  • Bernhard Grünewaldt
  • Elliotte Rusty Harold
  • Dennis Hendriks
  • Reid Hochstedler
  • Matthias Hummel
  • Curtis Jensen
  • Bruce A Johnson
  • Ismael Juma
  • Eugene Kirpichov
  • Oleksandr Kornieiev
  • Piotr Kochanski
  • Sergei Lebedev
  • Bob MacCallum
  • Jake Mannix
  • Benjamin McCann
  • Patrick Meyer
  • J. Lewis Muir
  • Venkatesha Murthy
  • Christopher Nix
  • Fredrik Norin
  • Sean Owen
  • Sujit Pal
  • Todd C. Parnell
  • Andreas Rieger
  • Sébastien Riou
  • Bill Rossi
  • Matthew Rowles
  • Pavel Ryzhov
  • Joni Salonen
  • Michael Saunders
  • Thorsten Schaefer
  • Christopher Schuck
  • Christian Semrau
  • David Stefka
  • Mauro Talevi
  • Radoslav Tsvetkov
  • Kim van der Linde
  • Alexey Volkov
  • Andrew Waterman
  • Jörg Weimar
  • Christian Winter
  • Piotr Wydrych
  • Xiaogang Zhang
  • Chris Popp
  • Mikkel Meyer Andersen
  • Bill Barker
  • Sébastien Brisard
  • Albert Davidson Chou
  • Mark Diggory
  • Robert Burrell Donkin
  • Otmar Ertl
  • Luc Maisonobe
  • Tim O'Brien
  • J. Pietschmann
  • Dimitri Pourbaix
  • Gilles Sadowski
  • Greg Sterijevski
  • Brent Worden
  • Thomas Neidhart
  • Evan Ward 
/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.math3.transform;

import java.io.Serializable;

import org.apache.commons.math3.analysis.FunctionUtils;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.ArithmeticUtils;


Implements the Fast Hadamard Transform (FHT).
Transformation of an input vector x to the output vector y.


In addition to transformation of real vectors, the Hadamard transform can
transform integer vectors into integer vectors. However, this integer transform
cannot be inverted directly. Due to a scaling factor it may lead to rational results.
As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
vector (1/2, -1/2, 0, 0).

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

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

    
{@inheritDoc}

Throws:
/**
     * {@inheritDoc}
     *
     * @throws MathIllegalArgumentException if the length of the data array is
     * not a power of two
     */
    public double[] transform(final double[] f, final TransformType type) {
        if (type == TransformType.FORWARD) {
            return fht(f);
        }
        return TransformUtils.scaleArray(fht(f), 1.0 / f.length);
    }

    
{@inheritDoc}

Throws:
/**
     * {@inheritDoc}
     *
     * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
     *   if the lower bound is greater than, or equal to the upper bound
     * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
     *   if the number of sample points is negative
     * @throws MathIllegalArgumentException if the number of sample points is not a power of two
     */
    public double[] transform(final UnivariateFunction f,
        final double min, final double max, final int n,
        final TransformType type) {

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

    
Returns the forward transform of the specified integer data set.The
integer transform cannot be inverted directly, due to a scaling factor
which may lead to double results.

Params:
  • f – the integer data array to be transformed (signal)
Throws:
Returns:the integer transformed array (spectrum)
/**
     * Returns the forward transform of the specified integer data set.The
     * integer transform cannot be inverted directly, due to a scaling factor
     * which may lead to double results.
     *
     * @param f the integer data array to be transformed (signal)
     * @return the integer transformed array (spectrum)
     * @throws MathIllegalArgumentException if the length of the data array is not a power of two
     */
    public int[] transform(final int[] f) {
        return fht(f);
    }

    
The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. Requires N * log2(N) = n * 2^n additions. 
Short Table of manual calculation for N=8


    

  1. x is the input vector to be transformed,
    2. 

  2. y is the output vector (Fast Hadamard transform of x),
    3. 

  3. a and b are helper rows.
    4. 







    x
    a
    b
    y




    x0
    a0 = x0 + x1
    b0 = a0 + a1
    y0 = b0+ b1




    x1
    a1 = x2 + x3
    b0 = a2 + a3
    y0 = b2 + b3




    x2
    a2 = x4 + x5
    b0 = a4 + a5
    y0 = b4 + b5




    x3
    a3 = x6 + x7
    b0 = a6 + a7
    y0 = b6 + b7




    x4
    a0 = x0 - x1
    b0 = a0 - a1
    y0 = b0 - b1




    x5
    a1 = x2 - x3
    b0 = a2 - a3
    y0 = b2 - b3




    x6
    a2 = x4 - x5
    b0 = a4 - a5
    y0 = b4 - b5




    x7
    a3 = x6 - x7
    b0 = a6 - a7
    y0 = b6 - b7






How it works


    

  1. Construct a matrix with N rows and n + 1 columns, hadm[n+1][N].
    
(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][N])
    2. 

  2. Place the input vector x[N] in the first column of the matrix hadm.
    3. 

  3. The entries of the submatrix D_top are calculated as follows 
      
        
    • D_top goes from entry [0][1] to [N / 2 - 1][n + 1],
      • 
        
    • the columns of D_top are the pairwise mutually exclusive sums of the previous column.
      • 
    
    

    4. 

  4. The entries of the submatrix D_bottom are calculated as follows 
      
        
    • D_bottom goes from entry [N / 2][1] to [N][n + 1],
      • 
        
    • the columns of D_bottom are the pairwise differences of the previous column.
      • 
    
    

    5. 

  5. The consputation of D_top and D_bottom are best understood with the above example (for N = 8). 
  6. The output vector y is now in the last column of hadm.
    7. 

  7. Algorithm from chipcenter.
    8. 



Visually






    0123
    …
    n + 1




    0
    x0
    
        ↑

        ← Dtop →

        ↓
    		



1x1


2x2


……


N / 2 - 1xN/2-1



    N / 2
    xN/2
    
        ↑

        ← Dbottom →

        ↓
    		



N / 2 + 1xN/2+1


N / 2 + 2xN/2+2


……


NxN





Params:
  • x – the real data array to be transformed
Throws:
Returns:the real transformed array, y
/**
     * The FHT (Fast Hadamard Transformation) which uses only subtraction and
     * addition. Requires {@code N * log2(N) = n * 2^n} additions.
     *
     * <h3>Short Table of manual calculation for N=8</h3>
     * <ol>
     * <li><b>x</b> is the input vector to be transformed,</li>
     * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>),</li>
     * <li>a and b are helper rows.</li>
     * </ol>
     * <table align="center" border="1" cellpadding="3">
     * <tbody align="center">
     * <tr>
     *     <th>x</th>
     *     <th>a</th>
     *     <th>b</th>
     *     <th>y</th>
     * </tr>
     * <tr>
     *     <th>x<sub>0</sub></th>
     *     <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td>
     *     <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td>
     *     <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>1</sub></th>
     *     <td>a<sub>1</sub> = x<sub>2</sub> + x<sub>3</sub></td>
     *     <td>b<sub>0</sub> = a<sub>2</sub> + a<sub>3</sub></td>
     *     <td>y<sub>0</sub> = b<sub>2</sub> + b<sub>3</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>2</sub></th>
     *     <td>a<sub>2</sub> = x<sub>4</sub> + x<sub>5</sub></td>
     *     <td>b<sub>0</sub> = a<sub>4</sub> + a<sub>5</sub></td>
     *     <td>y<sub>0</sub> = b<sub>4</sub> + b<sub>5</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>3</sub></th>
     *     <td>a<sub>3</sub> = x<sub>6</sub> + x<sub>7</sub></td>
     *     <td>b<sub>0</sub> = a<sub>6</sub> + a<sub>7</sub></td>
     *     <td>y<sub>0</sub> = b<sub>6</sub> + b<sub>7</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>4</sub></th>
     *     <td>a<sub>0</sub> = x<sub>0</sub> - x<sub>1</sub></td>
     *     <td>b<sub>0</sub> = a<sub>0</sub> - a<sub>1</sub></td>
     *     <td>y<sub>0</sub> = b<sub>0</sub> - b<sub>1</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>5</sub></th>
     *     <td>a<sub>1</sub> = x<sub>2</sub> - x<sub>3</sub></td>
     *     <td>b<sub>0</sub> = a<sub>2</sub> - a<sub>3</sub></td>
     *     <td>y<sub>0</sub> = b<sub>2</sub> - b<sub>3</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>6</sub></th>
     *     <td>a<sub>2</sub> = x<sub>4</sub> - x<sub>5</sub></td>
     *     <td>b<sub>0</sub> = a<sub>4</sub> - a<sub>5</sub></td>
     *     <td>y<sub>0</sub> = b<sub>4</sub> - b<sub>5</sub></td>
     * </tr>
     * <tr>
     *     <th>x<sub>7</sub></th>
     *     <td>a<sub>3</sub> = x<sub>6</sub> - x<sub>7</sub></td>
     *     <td>b<sub>0</sub> = a<sub>6</sub> - a<sub>7</sub></td>
     *     <td>y<sub>0</sub> = b<sub>6</sub> - b<sub>7</sub></td>
     * </tr>
     * </tbody>
     * </table>
     *
     * <h3>How it works</h3>
     * <ol>
     * <li>Construct a matrix with {@code N} rows and {@code n + 1} columns,
     * {@code hadm[n+1][N]}.<br/>
     * <em>(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][N])</em></li>
     * <li>Place the input vector {@code x[N]} in the first column of the
     * matrix {@code hadm}.</li>
     * <li>The entries of the submatrix {@code D_top} are calculated as follows
     *     <ul>
     *         <li>{@code D_top} goes from entry {@code [0][1]} to
     *         {@code [N / 2 - 1][n + 1]},</li>
     *         <li>the columns of {@code D_top} are the pairwise mutually
     *         exclusive sums of the previous column.</li>
     *     </ul>
     * </li>
     * <li>The entries of the submatrix {@code D_bottom} are calculated as
     * follows
     *     <ul>
     *         <li>{@code D_bottom} goes from entry {@code [N / 2][1]} to
     *         {@code [N][n + 1]},</li>
     *         <li>the columns of {@code D_bottom} are the pairwise differences
     *         of the previous column.</li>
     *     </ul>
     * </li>
     * <li>The consputation of {@code D_top} and {@code D_bottom} are best
     * understood with the above example (for {@code N = 8}).
     * <li>The output vector {@code y} is now in the last column of
     * {@code hadm}.</li>
     * <li><em>Algorithm from <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">chipcenter</a>.</em></li>
     * </ol>
     * <h3>Visually</h3>
     * <table border="1" align="center" cellpadding="3">
     * <tbody align="center">
     * <tr>
     *     <td></td><th>0</th><th>1</th><th>2</th><th>3</th>
     *     <th>&hellip;</th>
     *     <th>n + 1</th>
     * </tr>
     * <tr>
     *     <th>0</th>
     *     <td>x<sub>0</sub></td>
     *     <td colspan="5" rowspan="5" align="center" valign="middle">
     *         &uarr;<br/>
     *         &larr; D<sub>top</sub> &rarr;<br/>
     *         &darr;
     *     </td>
     * </tr>
     * <tr><th>1</th><td>x<sub>1</sub></td></tr>
     * <tr><th>2</th><td>x<sub>2</sub></td></tr>
     * <tr><th>&hellip;</th><td>&hellip;</td></tr>
     * <tr><th>N / 2 - 1</th><td>x<sub>N/2-1</sub></td></tr>
     * <tr>
     *     <th>N / 2</th>
     *     <td>x<sub>N/2</sub></td>
     *     <td colspan="5" rowspan="5" align="center" valign="middle">
     *         &uarr;<br/>
     *         &larr; D<sub>bottom</sub> &rarr;<br/>
     *         &darr;
     *     </td>
     * </tr>
     * <tr><th>N / 2 + 1</th><td>x<sub>N/2+1</sub></td></tr>
     * <tr><th>N / 2 + 2</th><td>x<sub>N/2+2</sub></td></tr>
     * <tr><th>&hellip;</th><td>&hellip;</td></tr>
     * <tr><th>N</th><td>x<sub>N</sub></td></tr>
     * </tbody>
     * </table>
     *
     * @param x the real data array to be transformed
     * @return the real transformed array, {@code y}
     * @throws MathIllegalArgumentException if the length of the data array is not a power of two
     */
    protected double[] fht(double[] x) throws MathIllegalArgumentException {

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

        if (!ArithmeticUtils.isPowerOfTwo(n)) {
            throw new MathIllegalArgumentException(
                    LocalizedFormats.NOT_POWER_OF_TWO,
                    Integer.valueOf(n));
        }

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

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

            // switch columns
            final double[] yTmp = yCurrent;
            yCurrent  = yPrevious;
            yPrevious = yTmp;

            // iterate from top to bottom (row)
            for (int i = 0; i < halfN; ++i) {
                // Dtop: the top part works with addition
                final int twoI = 2 * i;
                yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
            }
            for (int i = halfN; i < n; ++i) {
                // Dbottom: the bottom part works with subtraction
                final int twoI = 2 * i;
                yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
            }
        }

        return yCurrent;

    }

    
Returns the forward transform of the specified integer data set. The FHT
(Fast Hadamard Transform) uses only subtraction and addition.

Params:
  • x – the integer data array to be transformed
Throws:
Returns:the integer transformed array, y
/**
     * Returns the forward transform of the specified integer data set. The FHT
     * (Fast Hadamard Transform) uses only subtraction and addition.
     *
     * @param x the integer data array to be transformed
     * @return the integer transformed array, {@code y}
     * @throws MathIllegalArgumentException if the length of the data array is not a power of two
     */
    protected int[] fht(int[] x) throws MathIllegalArgumentException {

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

        if (!ArithmeticUtils.isPowerOfTwo(n)) {
            throw new MathIllegalArgumentException(
                    LocalizedFormats.NOT_POWER_OF_TWO,
                    Integer.valueOf(n));
        }

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

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

            // switch columns
            final int[] yTmp = yCurrent;
            yCurrent  = yPrevious;
            yPrevious = yTmp;

            // iterate from top to bottom (row)
            for (int i = 0; i < halfN; ++i) {
                // Dtop: the top part works with addition
                final int twoI = 2 * i;
                yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
            }
            for (int i = halfN; i < n; ++i) {
                // Dbottom: the bottom part works with subtraction
                final int twoI = 2 * i;
                yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
            }
        }

        // return the last computed output vector y
        return yCurrent;

    }

}