diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/transform/FastHadamardTransformer.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/transform/FastHadamardTransformer.java | 252 |
1 files changed, 252 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/transform/FastHadamardTransformer.java b/src/main/java/org/apache/commons/math/transform/FastHadamardTransformer.java new file mode 100644 index 0000000..db79c99 --- /dev/null +++ b/src/main/java/org/apache/commons/math/transform/FastHadamardTransformer.java @@ -0,0 +1,252 @@ +/* + * 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 <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).</p> + * @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. + * <p>The integer transform cannot be inverted directly, due to a scaling + * factor it may lead to double results.</p> + * @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. + * <br> + * Requires <b>Nlog2N = n2</b><sup>n</sup> additions. + * <br> + * <br> + * <b><u>Short Table of manual calculation for N=8:</u></b> + * <ol> + * <li><b>x</b> is the input vector we want to transform</li> + * <li><b>y</b> is the output vector which is our desired result</li> + * <li>a and b are just helper rows</li> + * </ol> + * <pre> + * <code> + * +----+----------+---------+----------+ + * | <b>x</b> | <b>a</b> | <b>b</b> | <b>y</b> | + * +----+----------+---------+----------+ + * | x<sub>0</sub> | a<sub>0</sub>=x<sub>0</sub>+x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>+a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>+b<sub>1</sub> | + * +----+----------+---------+----------+ + * | x<sub>1</sub> | a<sub>1</sub>=x<sub>2</sub>+x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>+a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>+b<sub>3</sub> | + * +----+----------+---------+----------+ + * | x<sub>2</sub> | a<sub>2</sub>=x<sub>4</sub>+x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>+a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>+b<sub>5</sub> | + * +----+----------+---------+----------+ + * | x<sub>3</sub> | a<sub>3</sub>=x<sub>6</sub>+x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>+a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>+b<sub>7</sub> | + * +----+----------+---------+----------+ + * | x<sub>4</sub> | a<sub>0</sub>=x<sub>0</sub>-x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>-a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>-b<sub>1</sub> | + * +----+----------+---------+----------+ + * | x<sub>5</sub> | a<sub>1</sub>=x<sub>2</sub>-x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>-a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>-b<sub>3</sub> | + * +----+----------+---------+----------+ + * | x<sub>6</sub> | a<sub>2</sub>=x<sub>4</sub>-x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>-a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>-b<sub>5</sub> | + * +----+----------+---------+----------+ + * | x<sub>7</sub> | a<sub>3</sub>=x<sub>6</sub>-x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>-a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>-b<sub>7</sub> | + * +----+----------+---------+----------+ + * </code> + * </pre> + * + * <b><u>How it works</u></b> + * <ol> + * <li>Construct a matrix with N rows and n+1 columns<br> <b>hadm[n+1][N]</b> + * <br><i>(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])</i></li> + * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li> + * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows. + * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1]. + * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column + * </li> + * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows. + * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1]. + * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column + * </li> + * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above. + * <li>The output vector y is now in the last column of <b>hadm</b></li> + * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li> + * </ol> + * <br> + * <b><u>Visually</u></b> + * <pre> + * +--------+---+---+---+-----+---+ + * | 0 | 1 | 2 | 3 | ... |n+1| + * +------+--------+---+---+---+-----+---+ + * |0 | x<sub>0</sub> | /\ | + * |1 | x<sub>1</sub> | || | + * |2 | x<sub>2</sub> | <= D<sub>top</sub> => | + * |... | ... | || | + * |N/2-1 | x<sub>N/2-1</sub> | \/ | + * +------+--------+---+---+---+-----+---+ + * |N/2 | x<sub>N/2</sub> | /\ | + * |N/2+1 | x<sub>N/2+1</sub> | || | + * |N/2+2 | x<sub>N/2+2</sub> | <= D<sub>bottom</sub> => | which is in the last column of the matrix + * |... | ... | || | + * |N | x<sub>N/2</sub> | \/ | + * +------+--------+---+---+---+-----+---+ + * </pre> + * + * @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) { + // D<sub>top</sub> + // 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) { + // D<sub>bottom</sub> + // 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) { + // D<sub>top</sub> + // 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) { + // D<sub>bottom</sub> + // 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; + + } + +} |