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 import org.apache.commons.math.exception.util.LocalizedFormats; 023 024 /** 025 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT). 026 * Transformation of an input vector x to the output vector y. 027 * <p>In addition to transformation of real vectors, the Hadamard transform can 028 * transform integer vectors into integer vectors. However, this integer transform 029 * cannot be inverted directly. Due to a scaling factor it may lead to rational results. 030 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational 031 * vector (1/2, -1/2, 0, 0).</p> 032 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $ 033 * @since 2.0 034 */ 035 public class FastHadamardTransformer implements RealTransformer { 036 037 /** {@inheritDoc} */ 038 public double[] transform(double f[]) 039 throws IllegalArgumentException { 040 return fht(f); 041 } 042 043 /** {@inheritDoc} */ 044 public double[] transform(UnivariateRealFunction f, 045 double min, double max, int n) 046 throws FunctionEvaluationException, IllegalArgumentException { 047 return fht(FastFourierTransformer.sample(f, min, max, n)); 048 } 049 050 /** {@inheritDoc} */ 051 public double[] inversetransform(double f[]) 052 throws IllegalArgumentException { 053 return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length); 054 } 055 056 /** {@inheritDoc} */ 057 public double[] inversetransform(UnivariateRealFunction f, 058 double min, double max, int n) 059 throws FunctionEvaluationException, IllegalArgumentException { 060 final double[] unscaled = 061 fht(FastFourierTransformer.sample(f, min, max, n)); 062 return FastFourierTransformer.scaleArray(unscaled, 1.0 / n); 063 } 064 065 /** 066 * Transform the given real data set. 067 * <p>The integer transform cannot be inverted directly, due to a scaling 068 * factor it may lead to double results.</p> 069 * @param f the integer data array to be transformed (signal) 070 * @return the integer transformed array (spectrum) 071 * @throws IllegalArgumentException if any parameters are invalid 072 */ 073 public int[] transform(int f[]) 074 throws IllegalArgumentException { 075 return fht(f); 076 } 077 078 /** 079 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition. 080 * <br> 081 * Requires <b>Nlog2N = n2</b><sup>n</sup> additions. 082 * <br> 083 * <br> 084 * <b><u>Short Table of manual calculation for N=8:</u></b> 085 * <ol> 086 * <li><b>x</b> is the input vector we want to transform</li> 087 * <li><b>y</b> is the output vector which is our desired result</li> 088 * <li>a and b are just helper rows</li> 089 * </ol> 090 * <pre> 091 * <code> 092 * +----+----------+---------+----------+ 093 * | <b>x</b> | <b>a</b> | <b>b</b> | <b>y</b> | 094 * +----+----------+---------+----------+ 095 * | 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> | 096 * +----+----------+---------+----------+ 097 * | 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> | 098 * +----+----------+---------+----------+ 099 * | 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 }