org.apache.commons.math.transform.FastHadamardTransformer Maven / Gradle / Ivy
/*
* 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.math.transform;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* 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).
* @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
* @since 2.0
*/
public class FastHadamardTransformer implements RealTransformer {
/** {@inheritDoc} */
public double[] transform(double f[])
throws IllegalArgumentException {
return fht(f);
}
/** {@inheritDoc} */
public double[] transform(UnivariateRealFunction f,
double min, double max, int n)
throws FunctionEvaluationException, IllegalArgumentException {
return fht(FastFourierTransformer.sample(f, min, max, n));
}
/** {@inheritDoc} */
public double[] inversetransform(double f[])
throws IllegalArgumentException {
return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
}
/** {@inheritDoc} */
public double[] inversetransform(UnivariateRealFunction f,
double min, double max, int n)
throws FunctionEvaluationException, IllegalArgumentException {
final double[] unscaled =
fht(FastFourierTransformer.sample(f, min, max, n));
return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
}
/**
* Transform the given real data set.
* The integer transform cannot be inverted directly, due to a scaling
* factor it may lead to double results.
* @param f the integer data array to be transformed (signal)
* @return the integer transformed array (spectrum)
* @throws IllegalArgumentException if any parameters are invalid
*/
public int[] transform(int f[])
throws IllegalArgumentException {
return fht(f);
}
/**
* The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
*
* Requires Nlog2N = n2n additions.
*
*
* Short Table of manual calculation for N=8:
*
* - x is the input vector we want to transform
* - y is the output vector which is our desired result
* - a and b are just 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 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][m])
* - Place the input vector x[N] in the first column of the matrix hadm
* - The entries of the submatrix Dtop are calculated as follows.
*
Dtop goes from entry [0][1] to [N/2-1][n+1].
*
The columns of Dtop are the pairwise mutually exclusive sums of the previous column
*
* - The entries of the submatrix Dbottom are calculated as follows.
*
Dbottom goes from entry [N/2][1] to [N][n+1].
*
The columns of Dbottom are the pairwise differences of the previous column
*
* - How Dtop and Dbottom you can understand best with the example for N=8 above.
*
- The output vector y is now in the last column of hadm
* - Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html
*
*
* Visually
*
* +--------+---+---+---+-----+---+
* | 0 | 1 | 2 | 3 | ... |n+1|
* +------+--------+---+---+---+-----+---+
* |0 | x0 | /\ |
* |1 | x1 | || |
* |2 | x2 | <= Dtop => |
* |... | ... | || |
* |N/2-1 | xN/2-1 | \/ |
* +------+--------+---+---+---+-----+---+
* |N/2 | xN/2 | /\ |
* |N/2+1 | xN/2+1 | || |
* |N/2+2 | xN/2+2 | <= Dbottom => | which is in the last column of the matrix
* |... | ... | || |
* |N | xN/2 | \/ |
* +------+--------+---+---+---+-----+---+
*
*
* @param x input vector
* @return y output vector
* @exception IllegalArgumentException if input array is not a power of 2
*/
protected double[] fht(double x[]) throws IllegalArgumentException {
// n is the row count of the input vector x
final int n = x.length;
final int halfN = n / 2;
// n has to be of the form n = 2^p !!
if (!FastFourierTransformer.isPowerOf2(n)) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.NOT_POWER_OF_TWO,
n);
}
// Instead of creating a matrix with p+1 columns and n rows
// we will use two single dimension arrays which we will use 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 the last computed output vector y
return yCurrent;
}
/**
* The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
* @param x input vector
* @return y output vector
* @exception IllegalArgumentException if input array is not a power of 2
*/
protected int[] fht(int x[]) throws IllegalArgumentException {
// n is the row count of the input vector x
final int n = x.length;
final int halfN = n / 2;
// n has to be of the form n = 2^p !!
if (!FastFourierTransformer.isPowerOf2(n)) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.NOT_POWER_OF_TWO,
n);
}
// Instead of creating a matrix with p+1 columns and n rows
// we will use two single dimension arrays which we will use 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;
}
}