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