001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math.transform; 018 019 import org.apache.commons.math.FunctionEvaluationException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.analysis.UnivariateRealFunction; 022 023 /** 024 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT). 025 * Transformation of an input vector x to the output vector y. 026 * <p>In addition to transformation of real vectors, the Hadamard transform can 027 * transform integer vectors into integer vectors. However, this integer transform 028 * cannot be inverted directly. Due to a scaling factor it may lead to rational results. 029 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational 030 * vector (1/2, -1/2, 0, 0).</p> 031 * @version $Revision: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 2009) $ 032 * @since 2.0 033 */ 034 public class FastHadamardTransformer implements RealTransformer { 035 036 /** {@inheritDoc} */ 037 public double[] transform(double f[]) 038 throws IllegalArgumentException { 039 return fht(f); 040 } 041 042 /** {@inheritDoc} */ 043 public double[] transform(UnivariateRealFunction f, 044 double min, double max, int n) 045 throws FunctionEvaluationException, IllegalArgumentException { 046 return fht(FastFourierTransformer.sample(f, min, max, n)); 047 } 048 049 /** {@inheritDoc} */ 050 public double[] inversetransform(double f[]) 051 throws IllegalArgumentException { 052 return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length); 053 } 054 055 /** {@inheritDoc} */ 056 public double[] inversetransform(UnivariateRealFunction f, 057 double min, double max, int n) 058 throws FunctionEvaluationException, IllegalArgumentException { 059 final double[] unscaled = 060 fht(FastFourierTransformer.sample(f, min, max, n)); 061 return FastFourierTransformer.scaleArray(unscaled, 1.0 / n); 062 } 063 064 /** 065 * Transform the given real data set. 066 * <p>The integer transform cannot be inverted directly, due to a scaling 067 * factor it may lead to double results.</p> 068 * @param f the integer data array to be transformed (signal) 069 * @return the integer transformed array (spectrum) 070 * @throws IllegalArgumentException if any parameters are invalid 071 */ 072 public int[] transform(int f[]) 073 throws IllegalArgumentException { 074 return fht(f); 075 } 076 077 /** 078 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. 079 * <br> 080 * Requires <b>Nlog2N = n2</b><sup>n</sup> additions. 081 * <br> 082 * <br> 083 * <b><u>Short Table of manual calculation for N=8:</u></b> 084 * <ol> 085 * <li><b>x</b> is the input vector we want to transform</li> 086 * <li><b>y</b> is the output vector which is our desired result</li> 087 * <li>a and b are just helper rows</li> 088 * </ol> 089 * <pre> 090 * <code> 091 * +----+----------+---------+----------+ 092 * | <b>x</b> | <b>a</b> | <b>b</b> | <b>y</b> | 093 * +----+----------+---------+----------+ 094 * | 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> | 095 * +----+----------+---------+----------+ 096 * | 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> | 097 * +----+----------+---------+----------+ 098 * | 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> | 099 * +----+----------+---------+----------+ 100 * | 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> | 101 * +----+----------+---------+----------+ 102 * | 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> | 103 * +----+----------+---------+----------+ 104 * | 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> | 105 * +----+----------+---------+----------+ 106 * | 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> | 107 * +----+----------+---------+----------+ 108 * | 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> | 109 * +----+----------+---------+----------+ 110 * </code> 111 * </pre> 112 * 113 * <b><u>How it works</u></b> 114 * <ol> 115 * <li>Construct a matrix with N rows and n+1 columns<br> <b>hadm[n+1][N]</b> 116 * <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> 117 * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li> 118 * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows. 119 * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1]. 120 * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column 121 * </li> 122 * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows. 123 * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1]. 124 * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column 125 * </li> 126 * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above. 127 * <li>The output vector y is now in the last column of <b>hadm</b></li> 128 * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li> 129 * </ol> 130 * <br> 131 * <b><u>Visually</u></b> 132 * <pre> 133 * +--------+---+---+---+-----+---+ 134 * | 0 | 1 | 2 | 3 | ... |n+1| 135 * +------+--------+---+---+---+-----+---+ 136 * |0 | x<sub>0</sub> | /\ | 137 * |1 | x<sub>1</sub> | || | 138 * |2 | x<sub>2</sub> | <= D<sub>top</sub> => | 139 * |... | ... | || | 140 * |N/2-1 | x<sub>N/2-1</sub> | \/ | 141 * +------+--------+---+---+---+-----+---+ 142 * |N/2 | x<sub>N/2</sub> | /\ | 143 * |N/2+1 | x<sub>N/2+1</sub> | || | 144 * |N/2+2 | x<sub>N/2+2</sub> | <= D<sub>bottom</sub> => | which is in the last column of the matrix 145 * |... | ... | || | 146 * |N | x<sub>N/2</sub> | \/ | 147 * +------+--------+---+---+---+-----+---+ 148 * </pre> 149 * 150 * @param x input vector 151 * @return y output vector 152 * @exception IllegalArgumentException if input array is not a power of 2 153 */ 154 protected double[] fht(double x[]) throws IllegalArgumentException { 155 156 // n is the row count of the input vector x 157 final int n = x.length; 158 final int halfN = n / 2; 159 160 // n has to be of the form n = 2^p !! 161 if (!FastFourierTransformer.isPowerOf2(n)) { 162 throw MathRuntimeException.createIllegalArgumentException( 163 "{0} is not a power of 2", 164 n); 165 } 166 167 // Instead of creating a matrix with p+1 columns and n rows 168 // we will use two single dimension arrays which we will use in an alternating way. 169 double[] yPrevious = new double[n]; 170 double[] yCurrent = x.clone(); 171 172 // iterate from left to right (column) 173 for (int j = 1; j < n; j <<= 1) { 174 175 // switch columns 176 final double[] yTmp = yCurrent; 177 yCurrent = yPrevious; 178 yPrevious = yTmp; 179 180 // iterate from top to bottom (row) 181 for (int i = 0; i < halfN; ++i) { 182 // D<sub>top</sub> 183 // The top part works with addition 184 final int twoI = 2 * i; 185 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; 186 } 187 for (int i = halfN; i < n; ++i) { 188 // D<sub>bottom</sub> 189 // The bottom part works with subtraction 190 final int twoI = 2 * i; 191 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; 192 } 193 } 194 195 // return the last computed output vector y 196 return yCurrent; 197 198 } 199 /** 200 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. 201 * @param x input vector 202 * @return y output vector 203 * @exception IllegalArgumentException if input array is not a power of 2 204 */ 205 protected int[] fht(int x[]) throws IllegalArgumentException { 206 207 // n is the row count of the input vector x 208 final int n = x.length; 209 final int halfN = n / 2; 210 211 // n has to be of the form n = 2^p !! 212 if (!FastFourierTransformer.isPowerOf2(n)) { 213 throw MathRuntimeException.createIllegalArgumentException( 214 "{0} is not a power of 2", 215 n); 216 } 217 218 // Instead of creating a matrix with p+1 columns and n rows 219 // we will use two single dimension arrays which we will use in an alternating way. 220 int[] yPrevious = new int[n]; 221 int[] yCurrent = x.clone(); 222 223 // iterate from left to right (column) 224 for (int j = 1; j < n; j <<= 1) { 225 226 // switch columns 227 final int[] yTmp = yCurrent; 228 yCurrent = yPrevious; 229 yPrevious = yTmp; 230 231 // iterate from top to bottom (row) 232 for (int i = 0; i < halfN; ++i) { 233 // D<sub>top</sub> 234 // The top part works with addition 235 final int twoI = 2 * i; 236 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; 237 } 238 for (int i = halfN; i < n; ++i) { 239 // D<sub>bottom</sub> 240 // The bottom part works with subtraction 241 final int twoI = 2 * i; 242 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; 243 } 244 } 245 246 // return the last computed output vector y 247 return yCurrent; 248 249 } 250 251 }