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