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.math4.transform;
18
19 import java.util.function.UnaryOperator;
20
21 import org.apache.commons.numbers.core.ArithmeticUtils;
22
23 /**
24 * <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
25 * <p>
26 * The FHT can also transform integer vectors into integer vectors.
27 * However, this transform cannot be inverted directly, due to a scaling
28 * factor that may lead to rational results (for example, the inverse
29 * transform of integer vector (0, 1, 0, 1) is vector (1/2, -1/2, 0, 0).
30 */
31 public class FastHadamardTransform implements RealTransform {
32 /** Operation to be performed. */
33 private final UnaryOperator<double[]> op;
34
35 /**
36 * Default constructor.
37 */
38 public FastHadamardTransform() {
39 this(false);
40 }
41
42 /**
43 * @param inverse Whether to perform the inverse transform.
44 */
45 public FastHadamardTransform(final boolean inverse) {
46 op = create(inverse);
47 }
48
49 /**
50 * {@inheritDoc}
51 *
52 * @throws IllegalArgumentException if the length of the data array is
53 * not a power of two.
54 */
55 @Override
56 public double[] apply(final double[] f) {
57 return op.apply(f);
58 }
59
60 /**
61 * Returns the forward transform of the given data set.
62 * The integer transform cannot be inverted directly, due to a scaling
63 * factor which may lead to double results.
64 *
65 * @param f Data array to be transformed (signal).
66 * @return the transformed array (spectrum).
67 * @throws IllegalArgumentException if the length of the data array is
68 * not a power of two.
69 */
70 public int[] apply(final int[] f) {
71 return fht(f);
72 }
73
74 /**
75 * FHT uses only subtraction and addition.
76 * It requires {@code N * log2(N) = n * 2^n} additions.
77 *
78 * <h3>Short Table of manual calculation for N=8</h3>
79 * <ol>
80 * <li><b>x</b> is the input vector to be transformed,</li>
81 * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>),</li>
82 * <li>a and b are helper rows.</li>
83 * </ol>
84 * <table style="text-align: center" border="">
85 * <caption>manual calculation for N=8</caption>
86 * <tbody style="text-align: center">
87 * <tr>
88 * <th>x</th>
89 * <th>a</th>
90 * <th>b</th>
91 * <th>y</th>
92 * </tr>
93 * <tr>
94 * <th>x<sub>0</sub></th>
95 * <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td>
96 * <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td>
97 * <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td>
98 * </tr>
99 * <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 }