[ 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.