/* * 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 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; import java.io.Serializable; /** * 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 */ public class FastHadamardTransformer implements RealTransformer, Serializable { /** Serializable version identifier. */ static final long serialVersionUID = 20120211L; /** * {@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 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. * * @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 * {@code N * log2(N) = n * 2^n} additions. * *

Short Table of manual calculation for N=8

* *

x is the input vector to be transformed, *
y is the output vector (Fast Hadamard transform of x), *
a and b are helper rows. *

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

x	a	b	y
x₀	a₀ = x₀ + x₁	b₀ = a₀ + a₁	y₀ = b₀+ b₁
x₁	a₁ = x₂ + x₃	b₀ = a₂ + a₃	y₀ = b₂ + b₃
x₂	a₂ = x₄ + x₅	b₀ = a₄ + a₅	y₀ = b₄ + b₅
x₃	a₃ = x₆ + x₇	b₀ = a₆ + a₇	y₀ = b₆ + b₇
x₄	a₀ = x₀ - x₁	b₀ = a₀ - a₁	y₀ = b₀ - b₁
x₅	a₁ = x₂ - x₃	b₀ = a₂ - a₃	y₀ = b₂ - b₃
x₆	a₂ = x₄ - x₅	b₀ = a₄ - a₅	y₀ = b₄ - b₅
x₇	a₃ = x₆ - x₇	b₀ = a₆ - a₇	y₀ = b₆ - b₇

* *

How it works

* *

Construct a matrix with {@code N} rows and {@code n + 1} columns, {@code 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]) *
Place the input vector {@code x[N]} in the first column of the matrix {@code hadm}. *
The entries of the submatrix {@code D_top} are calculated as follows *
- {@code D_top} goes from entry {@code [0][1]} to {@code [N / 2 - 1][n + 1]}, *
- the columns of {@code D_top} are the pairwise mutually exclusive sums of the * previous column. *
*
The entries of the submatrix {@code D_bottom} are calculated as follows *
- {@code D_bottom} goes from entry {@code [N / 2][1]} to {@code [N][n + 1]}, *
- the columns of {@code D_bottom} are the pairwise differences of the previous * column. *
*
The consputation of {@code D_top} and {@code D_bottom} are best understood with the * above example (for {@code N = 8}). *
The output vector {@code y} is now in the last column of {@code hadm}. *
Algorithm from chipcenter. *

* *

Visually

* * * * * * * * * * * * * * * * * * * * * * * * * * * *

	0	1	2	3	…	n + 1
0	x₀	* ↑ * ← D_top → * ↓ *
1	x₁
2	x₂
…	…
N / 2 - 1	x_N/2-1
N / 2	x_N/2	* ↑ * ← D_bottom → * ↓ *
N / 2 + 1	x_N/2+1
N / 2 + 2	x_N/2+2
…	…
N	x_N

* * @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. * * @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; } }