1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math.transform; 18 19 import org.apache.commons.math.FunctionEvaluationException; 20 import org.apache.commons.math.MathRuntimeException; 21 import org.apache.commons.math.analysis.UnivariateRealFunction; 22 import org.apache.commons.math.exception.util.LocalizedFormats; 23 24 /** 25 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT). 26 * Transformation of an input vector x to the output vector y. 27 * <p>In addition to transformation of real vectors, the Hadamard transform can 28 * transform integer vectors into integer vectors. However, this integer transform 29 * cannot be inverted directly. Due to a scaling factor it may lead to rational results. 30 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational 31 * vector (1/2, -1/2, 0, 0).</p> 32 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $ 33 * @since 2.0 34 */ 35 public class FastHadamardTransformer implements RealTransformer { 36 37 /** {@inheritDoc} */ 38 public double[] transform(double f[]) 39 throws IllegalArgumentException { 40 return fht(f); 41 } 42 43 /** {@inheritDoc} */ 44 public double[] transform(UnivariateRealFunction f, 45 double min, double max, int n) 46 throws FunctionEvaluationException, IllegalArgumentException { 47 return fht(FastFourierTransformer.sample(f, min, max, n)); 48 } 49 50 /** {@inheritDoc} */ 51 public double[] inversetransform(double f[]) 52 throws IllegalArgumentException { 53 return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length); 54 } 55 56 /** {@inheritDoc} */ 57 public double[] inversetransform(UnivariateRealFunction f, 58 double min, double max, int n) 59 throws FunctionEvaluationException, IllegalArgumentException { 60 final double[] unscaled = 61 fht(FastFourierTransformer.sample(f, min, max, n)); 62 return FastFourierTransformer.scaleArray(unscaled, 1.0 / n); 63 } 64 65 /** 66 * Transform the given real data set. 67 * <p>The integer transform cannot be inverted directly, due to a scaling 68 * factor it may lead to double results.</p> 69 * @param f the integer data array to be transformed (signal) 70 * @return the integer transformed array (spectrum) 71 * @throws IllegalArgumentException if any parameters are invalid 72 */ 73 public int[] transform(int f[]) 74 throws IllegalArgumentException { 75 return fht(f); 76 } 77 78 /** 79 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. 80 * <br> 81 * Requires <b>Nlog2N = n2</b><sup>n</sup> additions. 82 * <br> 83 * <br> 84 * <b><u>Short Table of manual calculation for N=8:</u></b> 85 * <ol> 86 * <li><b>x</b> is the input vector we want to transform</li> 87 * <li><b>y</b> is the output vector which is our desired result</li> 88 * <li>a and b are just helper rows</li> 89 * </ol> 90 * <pre> 91 * <code> 92 * +----+----------+---------+----------+ 93 * | <b>x</b> | <b>a</b> | <b>b</b> | <b>y</b> | 94 * +----+----------+---------+----------+ 95 * | 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> | 96 * +----+----------+---------+----------+ 97 * | 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> | 98 * +----+----------+---------+----------+ 99 * | 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> | 100 * +----+----------+---------+----------+ 101 * | 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> | 102 * +----+----------+---------+----------+ 103 * | 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> | 104 * +----+----------+---------+----------+ 105 * | 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> | 106 * +----+----------+---------+----------+ 107 * | 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> | 108 * +----+----------+---------+----------+ 109 * | 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> | 110 * +----+----------+---------+----------+ 111 * </code> 112 * </pre> 113 * 114 * <b><u>How it works</u></b> 115 * <ol> 116 * <li>Construct a matrix with N rows and n+1 columns<br> <b>hadm[n+1][N]</b> 117 * <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> 118 * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li> 119 * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows. 120 * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1]. 121 * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column 122 * </li> 123 * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows. 124 * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1]. 125 * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column 126 * </li> 127 * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above. 128 * <li>The output vector y is now in the last column of <b>hadm</b></li> 129 * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li> 130 * </ol> 131 * <br> 132 * <b><u>Visually</u></b> 133 * <pre> 134 * +--------+---+---+---+-----+---+ 135 * | 0 | 1 | 2 | 3 | ... |n+1| 136 * +------+--------+---+---+---+-----+---+ 137 * |0 | x<sub>0</sub> | /\ | 138 * |1 | x<sub>1</sub> | || | 139 * |2 | x<sub>2</sub> | <= D<sub>top</sub> => | 140 * |... | ... | || | 141 * |N/2-1 | x<sub>N/2-1</sub> | \/ | 142 * +------+--------+---+---+---+-----+---+ 143 * |N/2 | x<sub>N/2</sub> | /\ | 144 * |N/2+1 | x<sub>N/2+1</sub> | || | 145 * |N/2+2 | x<sub>N/2+2</sub> | <= D<sub>bottom</sub> => | which is in the last column of the matrix 146 * |... | ... | || | 147 * |N | x<sub>N/2</sub> | \/ | 148 * +------+--------+---+---+---+-----+---+ 149 * </pre> 150 * 151 * @param x input vector 152 * @return y output vector 153 * @exception IllegalArgumentException if input array is not a power of 2 154 */ 155 protected double[] fht(double x[]) throws IllegalArgumentException { 156 157 // n is the row count of the input vector x 158 final int n = x.length; 159 final int halfN = n / 2; 160 161 // n has to be of the form n = 2^p !! 162 if (!FastFourierTransformer.isPowerOf2(n)) { 163 throw MathRuntimeException.createIllegalArgumentException( 164 LocalizedFormats.NOT_POWER_OF_TWO, 165 n); 166 } 167 168 // Instead of creating a matrix with p+1 columns and n rows 169 // we will use two single dimension arrays which we will use in an alternating way. 170 double[] yPrevious = new double[n]; 171 double[] yCurrent = x.clone(); 172 173 // iterate from left to right (column) 174 for (int j = 1; j < n; j <<= 1) { 175 176 // switch columns 177 final double[] yTmp = yCurrent; 178 yCurrent = yPrevious; 179 yPrevious = yTmp; 180 181 // iterate from top to bottom (row) 182 for (int i = 0; i < halfN; ++i) { 183 // D<sub>top</sub> 184 // The top part works with addition 185 final int twoI = 2 * i; 186 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; 187 } 188 for (int i = halfN; i < n; ++i) { 189 // D<sub>bottom</sub> 190 // The bottom part works with subtraction 191 final int twoI = 2 * i; 192 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; 193 } 194 } 195 196 // return the last computed output vector y 197 return yCurrent; 198 199 } 200 /** 201 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. 202 * @param x input vector 203 * @return y output vector 204 * @exception IllegalArgumentException if input array is not a power of 2 205 */ 206 protected int[] fht(int x[]) throws IllegalArgumentException { 207 208 // n is the row count of the input vector x 209 final int n = x.length; 210 final int halfN = n / 2; 211 212 // n has to be of the form n = 2^p !! 213 if (!FastFourierTransformer.isPowerOf2(n)) { 214 throw MathRuntimeException.createIllegalArgumentException( 215 LocalizedFormats.NOT_POWER_OF_TWO, 216 n); 217 } 218 219 // Instead of creating a matrix with p+1 columns and n rows 220 // we will use two single dimension arrays which we will use in an alternating way. 221 int[] yPrevious = new int[n]; 222 int[] yCurrent = x.clone(); 223 224 // iterate from left to right (column) 225 for (int j = 1; j < n; j <<= 1) { 226 227 // switch columns 228 final int[] yTmp = yCurrent; 229 yCurrent = yPrevious; 230 yPrevious = yTmp; 231 232 // iterate from top to bottom (row) 233 for (int i = 0; i < halfN; ++i) { 234 // D<sub>top</sub> 235 // The top part works with addition 236 final int twoI = 2 * i; 237 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; 238 } 239 for (int i = halfN; i < n; ++i) { 240 // D<sub>bottom</sub> 241 // The bottom part works with subtraction 242 final int twoI = 2 * i; 243 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; 244 } 245 } 246 247 // return the last computed output vector y 248 return yCurrent; 249 250 } 251 252 } 253
