You are viewing a plain text version of this content. The canonical link for it is here.
Posted to issues@commons.apache.org by "Bernhard Grünewaldt (JIRA)" <ji...@apache.org> on 2008/12/27 20:36:44 UTC
[jira] Issue Comment Edited: (MATH-215) Fast Hadamard Transform
[ https://issues.apache.org/jira/browse/MATH-215?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12659277#action_12659277 ]
gruenebe edited comment on MATH-215 at 12/27/08 11:35 AM:
---------------------------------------------------------------------
Ok, here is my tested first draft.
What needs to be done is test corner cases, but since I don't really know enough about this someone else should write the corner case tests.
A method for int vectors would be good, so that there are no floatingpoint-operations.
But before doing this I would at first like to hear your opinions.
{code:title=org.apache.commons.math.transform.FastHadamardTransformer|borderStyle=solid}
/*
* 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 java.io.Serializable;
/**
* 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.
*
*/
public class FastHadamardTransformer implements Serializable {
private static final long serialVersionUID = 5044269102877526860L;
/**
* Wrapper method for fht() for double vectors
*
* @param x input vector
* @return y output vector
* @throws IllegalArgumentException
*/
public double[] transform(double x[]) throws IllegalArgumentException {
return fht(x);
}
/**
* 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
* @throws IllegalArgumentException
*/
protected double[] fht(double x[]) throws IllegalArgumentException {
// N is the row count of the input vector x
int N = x.length;
// Instead of creating a matrix with n+1 columns and N rows
// we will use two single dimension arrays which we will use in an alternating way.
// The method parameter x will be our first single dimension array, so we need just one more which is y.
// y will hold the final result.
double[] y = new double[N];
// calculate n via <i>n=log<sub>2</sub>(N)</i>
int n = (int) (Math.log( N ) / Math.log( 2 ));
// N has to be of the form N=2^n !!
if (N != Integer.highestOneBit(N)) {
throw new IllegalArgumentException("N has to be power of 2!");
}
// The algorithm description says: "place the input vector x in the first column"
// We are using x and y, and we will start with x as the first and y as the second vector
// so x is already the "first column"
// iterate from left to right (column)
for (int j=1;j<n+1;j++) {
// iterate from top to bottom (row)
for (int i=0;i<N;i++) {
if (i<N/2) {
// D<sub>top</sub>
// The top part works with addition
if (j % 2 == 1) {
// with j mod 2 = 1 -> x holds the result of the previous iteration
y[i] = x[i*2] + x[i*2 +1];
} else {
// with j mod 2 = 0 -> y holds the result of the previous iteration
x[i] = y[i*2] + y[i*2 +1];
}
} else {
// D<sub>bottom</sub>
// The bottom part works with subtraction
if (j % 2 == 1) {
// with j mod 2 = 1 -> x holds the result of the previous iteration
y[i] = x[(i-N/2)*2] - x[(i-N/2)*2 +1];
} else {
// with j mod 2 = 0 -> y holds the result of the previous iteration
x[i] = y[(i-N/2)*2] - y[(i-N/2)*2 +1];
}
}
}
}
// return the computed output vector y
return y;
}
}
{code}
{code:title=org.apache.commons.math.transform.FastHadamardTransformerTest|borderStyle=solid}
/*
* 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 junit.framework.TestCase;
/**
* JUnit Test for HadamardTransformerTest
* @see org.apache.commons.math.transform.FastHadamardTransformer
*/
public final class FastHadamardTransformerTest extends TestCase {
/**
* Test of transformer for the a 8-point FHT (means N=8)
*/
public void test1() {
// Initiate the transformer
FastHadamardTransformer transformer = new FastHadamardTransformer();
// input vector x
double x[] = {1.0,4.0,-2.0,3.0,0,1.0,4.0,-1.0};
// output vector y
double y[] = {10.0,-4.0,2.0,-4.0,2.0,-12.0,6.0,8.0};
// transform input vector x to output vector
double result[] = transformer.transform(x);
for (int i=0;i<result.length;i++) {
// compare computed results to precomputed results
assertEquals(y[i], result[i]);
}
}
}
{code}
was (Author: gruenebe):
Ok, here is my tested first draft.
What needs to be done is test corner cases, but since I don't really know enough about this someone else should write the corner case tests.
A method for int vectors would be good, so that there are no floatingpoint-operations.
But before doing this I would at first like to hear your opinions.
{code:title=org.apache.commons.math.transform.FastHadamardTransformer|borderStyle=solid}
/*
* 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 java.io.Serializable;
import org.apache.commons.math.MathException;
import org.apache.commons.math.linear.RealMatrixImpl;
/**
* Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
*/
public class FastHadamardTransformer implements Serializable {
private static final long serialVersionUID = 5044269102877526860L;
/**
* Wrapper method for fht() for double vectors
*
* @param x input vector
* @return y output vector
* @throws MathException
* @throws IllegalArgumentException
*/
public double[] transform(double x[]) throws IllegalArgumentException {
return fht(x);
}
/**
* 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> => |
* |... | ... | || |
* |N | x<sub>N/2</sub> | \/ |
* +------+--------+---+---+---+-----+---+
* </pre>
*
* @param x input vector
* @return y output vector
* @throws MathException
* @throws IllegalArgumentException
*/
protected double[] fht(double x[]) throws IllegalArgumentException {
// N is the row count of the input vector x
int N = x.length;
// calculate n via <i>n=log<sub>2</sub>(N)</i>
int n = (int) (Math.log( N ) / Math.log( 2 ));
// N has to be of the form N=2^n !!
if (N != Math.pow(2, n)) {
throw new IllegalArgumentException("The row count (array length) of the input vexctor x hast to be N=2^n with n is int!");
}
// create a matrix with n+1 colums and N rows
RealMatrixImpl hadm = new RealMatrixImpl(N,n+1);
// place the input vector x in the first column
hadm.setColumn(0, x);
// iterate over the matrix
for (int j=1;j<n+1;j++) {
for (int i=0;i<N;i++) {
if (i<N/2) {
// D<sub>top</sub>
// The top part works with addition
hadm.setEntry(i, j, hadm.getEntry(i*2, j-1) + hadm.getEntry(i*2 +1, j-1) );
} else {
// D<sub>bottom</sub>
// The bottom part works with subtraction
hadm.setEntry(i, j, hadm.getEntry((i-N/2)*2, j-1) - hadm.getEntry((i-N/2)*2 +1, j-1) );
}
}
}
// return the computed output vector y which is in the last column of the matrix
return hadm.getColumn(n);
}
}
{code}
{code:title=org.apache.commons.math.transform.FastHadamardTransformerTest|borderStyle=solid}
/*
* 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 junit.framework.TestCase;
import org.apache.commons.math.MathException;
/**
* JUnit Test for HadamardTransformerTest
* @see org.apache.commons.math.transform.FastHadamardTransformer
*/
public final class FastHadamardTransformerTest extends TestCase {
/**
* Test of transformer for the a 8-point FHT (means N=8)
*/
public void test1() throws MathException {
// Initiate the transformer
FastHadamardTransformer transformer = new FastHadamardTransformer();
// input vector x
double x[] = {1.0,4.0,-2.0,3.0,0,1.0,4.0,-1.0};
// output vector y
double y[] = {10.0,-4.0,2.0,-4.0,2.0,-12.0,6.0,8.0};
// transform input vector x to output vector
double result[] = transformer.transform(x);
for (int i=0;i<result.length;i++) {
// compare computed results to precomputed results
assertEquals(y[i], result[i]);
}
}
}
{code}
> Fast Hadamard Transform
> -----------------------
>
> Key: MATH-215
> URL: https://issues.apache.org/jira/browse/MATH-215
> Project: Commons Math
> Issue Type: New Feature
> Affects Versions: 1.0, 1.1, 1.2
> Reporter: Daniel Kuan
> Fix For: 2.1
>
> Attachments: FastHadamardTransformer.java.diff, FastHadamardTransformerTest.java.diff
>
>
> To date, the mathematical transforms package of Commons Maths, org.apache.commons.math.transform, only contains implementations for the Fourier, Sine, and Cosine transforms.
> This issue serves to propose and track the creation of an implementation for the Hadamard transform.
> Definition of the hadamard transform:
> http://en.wikipedia.org/wiki/Hadamard_transform#Definition
> Unfortunately, Mathworld does not provide a very detailed definition.
> http://mathworld.wolfram.com/HadamardTransform.html
> An elegant algorithm for the fast hadamard transform can be found here:
> http://www.archive.chipcenter.com/dsp/DSP000517F1.html
--
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.