Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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 févr. 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 | } |