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 */ 017package org.apache.commons.math3.transform; 018 019import java.io.Serializable; 020 021import org.apache.commons.math3.analysis.FunctionUtils; 022import org.apache.commons.math3.analysis.UnivariateFunction; 023import org.apache.commons.math3.exception.MathIllegalArgumentException; 024import org.apache.commons.math3.exception.util.LocalizedFormats; 025import org.apache.commons.math3.util.ArithmeticUtils; 026 027/** 028 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT). 029 * Transformation of an input vector x to the output vector y. 030 * <p> 031 * In addition to transformation of real vectors, the Hadamard transform can 032 * transform integer vectors into integer vectors. However, this integer transform 033 * cannot be inverted directly. Due to a scaling factor it may lead to rational results. 034 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational 035 * vector (1/2, -1/2, 0, 0). 036 * 037 * @since 2.0 038 */ 039public class FastHadamardTransformer implements RealTransformer, Serializable { 040 041 /** Serializable version identifier. */ 042 static final long serialVersionUID = 20120211L; 043 044 /** 045 * {@inheritDoc} 046 * 047 * @throws MathIllegalArgumentException if the length of the data array is 048 * not a power of two 049 */ 050 public double[] transform(final double[] f, final TransformType type) { 051 if (type == TransformType.FORWARD) { 052 return fht(f); 053 } 054 return TransformUtils.scaleArray(fht(f), 1.0 / f.length); 055 } 056 057 /** 058 * {@inheritDoc} 059 * 060 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException 061 * if the lower bound is greater than, or equal to the upper bound 062 * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException 063 * if the number of sample points is negative 064 * @throws MathIllegalArgumentException if the number of sample points is not a power of two 065 */ 066 public double[] transform(final UnivariateFunction f, 067 final double min, final double max, final int n, 068 final TransformType type) { 069 070 return transform(FunctionUtils.sample(f, min, max, n), type); 071 } 072 073 /** 074 * Returns the forward transform of the specified integer data set.The 075 * integer transform cannot be inverted directly, due to a scaling factor 076 * which may lead to double results. 077 * 078 * @param f the integer data array to be transformed (signal) 079 * @return the integer transformed array (spectrum) 080 * @throws MathIllegalArgumentException if the length of the data array is not a power of two 081 */ 082 public int[] transform(final int[] f) { 083 return fht(f); 084 } 085 086 /** 087 * The FHT (Fast Hadamard Transformation) which uses only subtraction and 088 * addition. Requires {@code N * log2(N) = n * 2^n} additions. 089 * 090 * <h3>Short Table of manual calculation for N=8</h3> 091 * <ol> 092 * <li><b>x</b> is the input vector to be transformed,</li> 093 * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>),</li> 094 * <li>a and b are helper rows.</li> 095 * </ol> 096 * <table align="center" border="1" cellpadding="3"> 097 * <tbody align="center"> 098 * <tr> 099 * <th>x</th> 100 * <th>a</th> 101 * <th>b</th> 102 * <th>y</th> 103 * </tr> 104 * <tr> 105 * <th>x<sub>0</sub></th> 106 * <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td> 107 * <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td> 108 * <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td> 109 * </tr> 110 * <tr> 111 * <th>x<sub>1</sub></th> 112 * <td>a<sub>1</sub> = x<sub>2</sub> + x<sub>3</sub></td> 113 * <td>b<sub>0</sub> = a<sub>2</sub> + a<sub>3</sub></td> 114 * <td>y<sub>0</sub> = b<sub>2</sub> + b<sub>3</sub></td> 115 * </tr> 116 * <tr> 117 * <th>x<sub>2</sub></th> 118 * <td>a<sub>2</sub> = x<sub>4</sub> + x<sub>5</sub></td> 119 * <td>b<sub>0</sub> = a<sub>4</sub> + a<sub>5</sub></td> 120 * <td>y<sub>0</sub> = b<sub>4</sub> + b<sub>5</sub></td> 121 * </tr> 122 * <tr> 123 * <th>x<sub>3</sub></th> 124 * <td>a<sub>3</sub> = x<sub>6</sub> + x<sub>7</sub></td> 125 * <td>b<sub>0</sub> = a<sub>6</sub> + a<sub>7</sub></td> 126 * <td>y<sub>0</sub> = b<sub>6</sub> + b<sub>7</sub></td> 127 * </tr> 128 * <tr> 129 * <th>x<sub>4</sub></th> 130 * <td>a<sub>0</sub> = x<sub>0</sub> - x<sub>1</sub></td> 131 * <td>b<sub>0</sub> = a<sub>0</sub> - a<sub>1</sub></td> 132 * <td>y<sub>0</sub> = b<sub>0</sub> - b<sub>1</sub></td> 133 * </tr> 134 * <tr> 135 * <th>x<sub>5</sub></th> 136 * <td>a<sub>1</sub> = x<sub>2</sub> - x<sub>3</sub></td> 137 * <td>b<sub>0</sub> = a<sub>2</sub> - a<sub>3</sub></td> 138 * <td>y<sub>0</sub> = b<sub>2</sub> - b<sub>3</sub></td> 139 * </tr> 140 * <tr> 141 * <th>x<sub>6</sub></th> 142 * <td>a<sub>2</sub> = x<sub>4</sub> - x<sub>5</sub></td> 143 * <td>b<sub>0</sub> = a<sub>4</sub> - a<sub>5</sub></td> 144 * <td>y<sub>0</sub> = b<sub>4</sub> - b<sub>5</sub></td> 145 * </tr> 146 * <tr> 147 * <th>x<sub>7</sub></th> 148 * <td>a<sub>3</sub> = x<sub>6</sub> - x<sub>7</sub></td> 149 * <td>b<sub>0</sub> = a<sub>6</sub> - a<sub>7</sub></td> 150 * <td>y<sub>0</sub> = b<sub>6</sub> - b<sub>7</sub></td> 151 * </tr> 152 * </tbody> 153 * </table> 154 * 155 * <h3>How it works</h3> 156 * <ol> 157 * <li>Construct a matrix with {@code N} rows and {@code n + 1} columns, 158 * {@code hadm[n+1][N]}.<br/> 159 * <em>(If I use [x][y] it always means [row-offset][column-offset] of a 160 * Matrix with n rows and m columns. Its entries go from M[0][0] 161 * to M[n][N])</em></li> 162 * <li>Place the input vector {@code x[N]} in the first column of the 163 * matrix {@code hadm}.</li> 164 * <li>The entries of the submatrix {@code D_top} are calculated as follows 165 * <ul> 166 * <li>{@code D_top} goes from entry {@code [0][1]} to 167 * {@code [N / 2 - 1][n + 1]},</li> 168 * <li>the columns of {@code D_top} are the pairwise mutually 169 * exclusive sums of the previous column.</li> 170 * </ul> 171 * </li> 172 * <li>The entries of the submatrix {@code D_bottom} are calculated as 173 * follows 174 * <ul> 175 * <li>{@code D_bottom} goes from entry {@code [N / 2][1]} to 176 * {@code [N][n + 1]},</li> 177 * <li>the columns of {@code D_bottom} are the pairwise differences 178 * of the previous column.</li> 179 * </ul> 180 * </li> 181 * <li>The consputation of {@code D_top} and {@code D_bottom} are best 182 * understood with the above example (for {@code N = 8}). 183 * <li>The output vector {@code y} is now in the last column of 184 * {@code hadm}.</li> 185 * <li><em>Algorithm from <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">chipcenter</a>.</em></li> 186 * </ol> 187 * <h3>Visually</h3> 188 * <table border="1" align="center" cellpadding="3"> 189 * <tbody align="center"> 190 * <tr> 191 * <td></td><th>0</th><th>1</th><th>2</th><th>3</th> 192 * <th>…</th> 193 * <th>n + 1</th> 194 * </tr> 195 * <tr> 196 * <th>0</th> 197 * <td>x<sub>0</sub></td> 198 * <td colspan="5" rowspan="5" align="center" valign="middle"> 199 * ↑<br/> 200 * ← D<sub>top</sub> →<br/> 201 * ↓ 202 * </td> 203 * </tr> 204 * <tr><th>1</th><td>x<sub>1</sub></td></tr> 205 * <tr><th>2</th><td>x<sub>2</sub></td></tr> 206 * <tr><th>…</th><td>…</td></tr> 207 * <tr><th>N / 2 - 1</th><td>x<sub>N/2-1</sub></td></tr> 208 * <tr> 209 * <th>N / 2</th> 210 * <td>x<sub>N/2</sub></td> 211 * <td colspan="5" rowspan="5" align="center" valign="middle"> 212 * ↑<br/> 213 * ← D<sub>bottom</sub> →<br/> 214 * ↓ 215 * </td> 216 * </tr> 217 * <tr><th>N / 2 + 1</th><td>x<sub>N/2+1</sub></td></tr> 218 * <tr><th>N / 2 + 2</th><td>x<sub>N/2+2</sub></td></tr> 219 * <tr><th>…</th><td>…</td></tr> 220 * <tr><th>N</th><td>x<sub>N</sub></td></tr> 221 * </tbody> 222 * </table> 223 * 224 * @param x the real data array to be transformed 225 * @return the real transformed array, {@code y} 226 * @throws MathIllegalArgumentException if the length of the data array is not a power of two 227 */ 228 protected double[] fht(double[] x) throws MathIllegalArgumentException { 229 230 final int n = x.length; 231 final int halfN = n / 2; 232 233 if (!ArithmeticUtils.isPowerOfTwo(n)) { 234 throw new MathIllegalArgumentException( 235 LocalizedFormats.NOT_POWER_OF_TWO, 236 Integer.valueOf(n)); 237 } 238 239 /* 240 * Instead of creating a matrix with p+1 columns and n rows, we use two 241 * one dimension arrays which we are used in an alternating way. 242 */ 243 double[] yPrevious = new double[n]; 244 double[] yCurrent = x.clone(); 245 246 // iterate from left to right (column) 247 for (int j = 1; j < n; j <<= 1) { 248 249 // switch columns 250 final double[] yTmp = yCurrent; 251 yCurrent = yPrevious; 252 yPrevious = yTmp; 253 254 // iterate from top to bottom (row) 255 for (int i = 0; i < halfN; ++i) { 256 // Dtop: the top part works with addition 257 final int twoI = 2 * i; 258 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; 259 } 260 for (int i = halfN; i < n; ++i) { 261 // Dbottom: the bottom part works with subtraction 262 final int twoI = 2 * i; 263 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; 264 } 265 } 266 267 return yCurrent; 268 269 } 270 271 /** 272 * Returns the forward transform of the specified integer data set. The FHT 273 * (Fast Hadamard Transform) uses only subtraction and addition. 274 * 275 * @param x the integer data array to be transformed 276 * @return the integer transformed array, {@code y} 277 * @throws MathIllegalArgumentException if the length of the data array is not a power of two 278 */ 279 protected int[] fht(int[] x) throws MathIllegalArgumentException { 280 281 final int n = x.length; 282 final int halfN = n / 2; 283 284 if (!ArithmeticUtils.isPowerOfTwo(n)) { 285 throw new MathIllegalArgumentException( 286 LocalizedFormats.NOT_POWER_OF_TWO, 287 Integer.valueOf(n)); 288 } 289 290 /* 291 * Instead of creating a matrix with p+1 columns and n rows, we use two 292 * one dimension arrays which we are used in an alternating way. 293 */ 294 int[] yPrevious = new int[n]; 295 int[] yCurrent = x.clone(); 296 297 // iterate from left to right (column) 298 for (int j = 1; j < n; j <<= 1) { 299 300 // switch columns 301 final int[] yTmp = yCurrent; 302 yCurrent = yPrevious; 303 yPrevious = yTmp; 304 305 // iterate from top to bottom (row) 306 for (int i = 0; i < halfN; ++i) { 307 // Dtop: the top part works with addition 308 final int twoI = 2 * i; 309 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1]; 310 } 311 for (int i = halfN; i < n; ++i) { 312 // Dbottom: the bottom part works with subtraction 313 final int twoI = 2 * i; 314 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1]; 315 } 316 } 317 318 // return the last computed output vector y 319 return yCurrent; 320 321 } 322 323}