/*
* 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.exception.MathUserException;
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 $Id: FastHadamardTransformer.java 1131229 2011-06-03 20:49:25Z luc $
* @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 MathUserException, 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 MathUserException, 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;
}
}