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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.transform;
23
24 import java.io.Serializable;
25
26 import org.hipparchus.analysis.FunctionUtils;
27 import org.hipparchus.analysis.UnivariateFunction;
28 import org.hipparchus.exception.MathIllegalArgumentException;
29 import org.hipparchus.util.ArithmeticUtils;
30
31 /**
32 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
33 * Transformation of an input vector x to the output vector y.
34 * <p>
35 * In addition to transformation of real vectors, the Hadamard transform can
36 * transform integer vectors into integer vectors. However, this integer transform
37 * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
38 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
39 * vector (1/2, -1/2, 0, 0).
40 *
41 */
42 public class FastHadamardTransformer implements RealTransformer, Serializable {
43
44 /** Serializable version identifier. */
45 static final long serialVersionUID = 20120211L;
46
47 /** Empty constructor.
48 * <p>
49 * This constructor is not strictly necessary, but it prevents spurious
50 * javadoc warnings with JDK 18 and later.
51 * </p>
52 * @since 3.0
53 */
54 public FastHadamardTransformer() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
55 // nothing to do
56 }
57
58 /**
59 * {@inheritDoc}
60 *
61 * @throws MathIllegalArgumentException if the length of the data array is
62 * not a power of two
63 */
64 @Override
65 public double[] transform(final double[] f, final TransformType type) {
66 if (type == TransformType.FORWARD) {
67 return fht(f);
68 }
69 return TransformUtils.scaleArray(fht(f), 1.0 / f.length);
70 }
71
72 /**
73 * {@inheritDoc}
74 *
75 * @throws org.hipparchus.exception.MathIllegalArgumentException
76 * if the lower bound is greater than, or equal to the upper bound
77 * @throws org.hipparchus.exception.MathIllegalArgumentException
78 * if the number of sample points is negative
79 * @throws MathIllegalArgumentException if the number of sample points is not a power of two
80 */
81 @Override
82 public double[] transform(final UnivariateFunction f,
83 final double min, final double max, final int n,
84 final TransformType type) {
85
86 return transform(FunctionUtils.sample(f, min, max, n), type);
87 }
88
89 /**
90 * Returns the forward transform of the specified integer data set.The
91 * integer transform cannot be inverted directly, due to a scaling factor
92 * which may lead to double results.
93 *
94 * @param f the integer data array to be transformed (signal)
95 * @return the integer transformed array (spectrum)
96 * @throws MathIllegalArgumentException if the length of the data array is not a power of two
97 */
98 public int[] transform(final int[] f) {
99 return fht(f);
100 }
101
102 /**
103 * The FHT (Fast Hadamard Transformation) which uses only subtraction and
104 * addition. Requires {@code N * log2(N) = n * 2^n} additions.
105 *
106 * <ol>
107 * <li><b>x</b> is the input vector to be transformed,</li>
108 * <li><b>y</b> is the output vector (Fast Hadamard transform of <b>x</b>),</li>
109 * <li>a and b are helper rows.</li>
110 * </ol>
111 * <table border="1">
112 * <caption>Short Table of manual calculation for N=8</caption>
113 * <tbody>
114 * <tr>
115 * <th>x</th>
116 * <th>a</th>
117 * <th>b</th>
118 * <th>y</th>
119 * </tr>
120 * <tr>
121 * <th>x<sub>0</sub></th>
122 * <td>a<sub>0</sub> = x<sub>0</sub> + x<sub>1</sub></td>
123 * <td>b<sub>0</sub> = a<sub>0</sub> + a<sub>1</sub></td>
124 * <td>y<sub>0</sub> = b<sub>0</sub >+ b<sub>1</sub></td>
125 * </tr>
126 * <tr>
127 * <th>x<sub>1</sub></th>
128 * <td>a<sub>1</sub> = x<sub>2</sub> + x<sub>3</sub></td>
129 * <td>b<sub>0</sub> = a<sub>2</sub> + a<sub>3</sub></td>
130 * <td>y<sub>0</sub> = b<sub>2</sub> + b<sub>3</sub></td>
131 * </tr>
132 * <tr>
133 * <th>x<sub>2</sub></th>
134 * <td>a<sub>2</sub> = x<sub>4</sub> + x<sub>5</sub></td>
135 * <td>b<sub>0</sub> = a<sub>4</sub> + a<sub>5</sub></td>
136 * <td>y<sub>0</sub> = b<sub>4</sub> + b<sub>5</sub></td>
137 * </tr>
138 * <tr>
139 * <th>x<sub>3</sub></th>
140 * <td>a<sub>3</sub> = x<sub>6</sub> + x<sub>7</sub></td>
141 * <td>b<sub>0</sub> = a<sub>6</sub> + a<sub>7</sub></td>
142 * <td>y<sub>0</sub> = b<sub>6</sub> + b<sub>7</sub></td>
143 * </tr>
144 * <tr>
145 * <th>x<sub>4</sub></th>
146 * <td>a<sub>0</sub> = x<sub>0</sub> - x<sub>1</sub></td>
147 * <td>b<sub>0</sub> = a<sub>0</sub> - a<sub>1</sub></td>
148 * <td>y<sub>0</sub> = b<sub>0</sub> - b<sub>1</sub></td>
149 * </tr>
150 * <tr>
151 * <th>x<sub>5</sub></th>
152 * <td>a<sub>1</sub> = x<sub>2</sub> - x<sub>3</sub></td>
153 * <td>b<sub>0</sub> = a<sub>2</sub> - a<sub>3</sub></td>
154 * <td>y<sub>0</sub> = b<sub>2</sub> - b<sub>3</sub></td>
155 * </tr>
156 * <tr>
157 * <th>x<sub>6</sub></th>
158 * <td>a<sub>2</sub> = x<sub>4</sub> - x<sub>5</sub></td>
159 * <td>b<sub>0</sub> = a<sub>4</sub> - a<sub>5</sub></td>
160 * <td>y<sub>0</sub> = b<sub>4</sub> - b<sub>5</sub></td>
161 * </tr>
162 * <tr>
163 * <th>x<sub>7</sub></th>
164 * <td>a<sub>3</sub> = x<sub>6</sub> - x<sub>7</sub></td>
165 * <td>b<sub>0</sub> = a<sub>6</sub> - a<sub>7</sub></td>
166 * <td>y<sub>0</sub> = b<sub>6</sub> - b<sub>7</sub></td>
167 * </tr>
168 * </tbody>
169 * </table>
170 *
171 * <p>How it works</p>
172 * <ol>
173 * <li>Construct a matrix with {@code N} rows and {@code n + 1} columns,
174 * {@code hadm[n+1][N]}.<br>
175 * <em>(If I use [x][y] it always means [row-offset][column-offset] of a
176 * Matrix with n rows and m columns. Its entries go from M[0][0]
177 * to M[n][N])</em></li>
178 * <li>Place the input vector {@code x[N]} in the first column of the
179 * matrix {@code hadm}.</li>
180 * <li>The entries of the submatrix {@code D_top} are calculated as follows
181 * <ul>
182 * <li>{@code D_top} goes from entry {@code [0][1]} to
183 * {@code [N / 2 - 1][n + 1]},</li>
184 * <li>the columns of {@code D_top} are the pairwise mutually
185 * exclusive sums of the previous column.</li>
186 * </ul>
187 * </li>
188 * <li>The entries of the submatrix {@code D_bottom} are calculated as
189 * follows
190 * <ul>
191 * <li>{@code D_bottom} goes from entry {@code [N / 2][1]} to
192 * {@code [N][n + 1]},</li>
193 * <li>the columns of {@code D_bottom} are the pairwise differences
194 * of the previous column.</li>
195 * </ul>
196 * </li>
197 * <li>The consputation of {@code D_top} and {@code D_bottom} are best
198 * understood with the above example (for {@code N = 8}).
199 * <li>The output vector {@code y} is now in the last column of
200 * {@code hadm}.</li>
201 * <li><em>Algorithm from <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">chipcenter</a>.</em></li>
202 * </ol>
203 * <table border="1" >
204 * <caption>visually</caption>
205 * <tbody>
206 * <tr>
207 * <td></td><th>0</th><th>1</th><th>2</th><th>3</th>
208 * <th>…</th>
209 * <th>n + 1</th>
210 * </tr>
211 * <tr>
212 * <th>0</th>
213 * <td>x<sub>0</sub></td>
214 * <td colspan="5" rowspan="5" >
215 * ↑<br>
216 * ← D<sub>top</sub> →<br>
217 * ↓
218 * </td>
219 * </tr>
220 * <tr><th>1</th><td>x<sub>1</sub></td></tr>
221 * <tr><th>2</th><td>x<sub>2</sub></td></tr>
222 * <tr><th>…</th><td>…</td></tr>
223 * <tr><th>N / 2 - 1</th><td>x<sub>N/2-1</sub></td></tr>
224 * <tr>
225 * <th>N / 2</th>
226 * <td>x<sub>N/2</sub></td>
227 * <td colspan="5" rowspan="5" >
228 * ↑<br>
229 * ← D<sub>bottom</sub> →<br>
230 * ↓
231 * </td>
232 * </tr>
233 * <tr><th>N / 2 + 1</th><td>x<sub>N/2+1</sub></td></tr>
234 * <tr><th>N / 2 + 2</th><td>x<sub>N/2+2</sub></td></tr>
235 * <tr><th>…</th><td>…</td></tr>
236 * <tr><th>N</th><td>x<sub>N</sub></td></tr>
237 * </tbody>
238 * </table>
239 *
240 * @param x the real data array to be transformed
241 * @return the real transformed array, {@code y}
242 * @throws MathIllegalArgumentException if the length of the data array is not a power of two
243 */
244 protected double[] fht(double[] x) throws MathIllegalArgumentException {
245
246 final int n = x.length;
247 final int halfN = n / 2;
248
249 if (!ArithmeticUtils.isPowerOfTwo(n)) {
250 throw new MathIllegalArgumentException(LocalizedFFTFormats.NOT_POWER_OF_TWO,
251 Integer.valueOf(n));
252 }
253
254 /*
255 * Instead of creating a matrix with p+1 columns and n rows, we use two
256 * one dimension arrays which we are used in an alternating way.
257 */
258 double[] yPrevious = new double[n];
259 double[] yCurrent = x.clone();
260
261 // iterate from left to right (column)
262 for (int j = 1; j < n; j <<= 1) {
263
264 // switch columns
265 final double[] yTmp = yCurrent;
266 yCurrent = yPrevious;
267 yPrevious = yTmp;
268
269 // iterate from top to bottom (row)
270 for (int i = 0; i < halfN; ++i) {
271 // Dtop: the top part works with addition
272 final int twoI = 2 * i;
273 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
274 }
275 for (int i = halfN; i < n; ++i) {
276 // Dbottom: the bottom part works with subtraction
277 final int twoI = 2 * i;
278 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
279 }
280 }
281
282 return yCurrent;
283
284 }
285
286 /**
287 * Returns the forward transform of the specified integer data set. The FHT
288 * (Fast Hadamard Transform) uses only subtraction and addition.
289 *
290 * @param x the integer data array to be transformed
291 * @return the integer transformed array, {@code y}
292 * @throws MathIllegalArgumentException if the length of the data array is not a power of two
293 */
294 protected int[] fht(int[] x) throws MathIllegalArgumentException {
295
296 final int n = x.length;
297 final int halfN = n / 2;
298
299 if (!ArithmeticUtils.isPowerOfTwo(n)) {
300 throw new MathIllegalArgumentException(LocalizedFFTFormats.NOT_POWER_OF_TWO,
301 Integer.valueOf(n));
302 }
303
304 /*
305 * Instead of creating a matrix with p+1 columns and n rows, we use two
306 * one dimension arrays which we are used in an alternating way.
307 */
308 int[] yPrevious = new int[n];
309 int[] yCurrent = x.clone();
310
311 // iterate from left to right (column)
312 for (int j = 1; j < n; j <<= 1) {
313
314 // switch columns
315 final int[] yTmp = yCurrent;
316 yCurrent = yPrevious;
317 yPrevious = yTmp;
318
319 // iterate from top to bottom (row)
320 for (int i = 0; i < halfN; ++i) {
321 // Dtop: the top part works with addition
322 final int twoI = 2 * i;
323 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
324 }
325 for (int i = halfN; i < n; ++i) {
326 // Dbottom: the bottom part works with subtraction
327 final int twoI = 2 * i;
328 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
329 }
330 }
331
332 // return the last computed output vector y
333 return yCurrent;
334
335 }
336
337 }