org.apache.commons.math3.transform.FastHadamardTransformer Maven / Gradle / Ivy
Show all versions of virtdata-lib-realer Show documentation
/*
* 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
*/
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
*
*
* 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
*
* - 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
* x0
*
* ↑
* ← Dtop →
* ↓
*
*
* 1 x1 2 x2 … … N / 2 - 1 xN/2-1
* N / 2
* xN/2
*
* ↑
* ← Dbottom →
* ↓
*
*
* N / 2 + 1 xN/2+1 N / 2 + 2 xN/2+2 … … N xN
*
* @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;
}
}