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.math.transform;
18
19 import org.apache.commons.math.FunctionEvaluationException;
20 import org.apache.commons.math.MathRuntimeException;
21 import org.apache.commons.math.analysis.UnivariateRealFunction;
22
23 /**
24 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
25 * Transformation of an input vector x to the output vector y.
26 * <p>In addition to transformation of real vectors, the Hadamard transform can
27 * transform integer vectors into integer vectors. However, this integer transform
28 * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
29 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
30 * vector (1/2, -1/2, 0, 0).</p>
31 * @version $Revision: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 2009) $
32 * @since 2.0
33 */
34 public class FastHadamardTransformer implements RealTransformer {
35
36 /** {@inheritDoc} */
37 public double[] transform(double f[])
38 throws IllegalArgumentException {
39 return fht(f);
40 }
41
42 /** {@inheritDoc} */
43 public double[] transform(UnivariateRealFunction f,
44 double min, double max, int n)
45 throws FunctionEvaluationException, IllegalArgumentException {
46 return fht(FastFourierTransformer.sample(f, min, max, n));
47 }
48
49 /** {@inheritDoc} */
50 public double[] inversetransform(double f[])
51 throws IllegalArgumentException {
52 return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
53 }
54
55 /** {@inheritDoc} */
56 public double[] inversetransform(UnivariateRealFunction f,
57 double min, double max, int n)
58 throws FunctionEvaluationException, IllegalArgumentException {
59 final double[] unscaled =
60 fht(FastFourierTransformer.sample(f, min, max, n));
61 return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
62 }
63
64 /**
65 * Transform the given real data set.
66 * <p>The integer transform cannot be inverted directly, due to a scaling
67 * factor it may lead to double results.</p>
68 * @param f the integer data array to be transformed (signal)
69 * @return the integer transformed array (spectrum)
70 * @throws IllegalArgumentException if any parameters are invalid
71 */
72 public int[] transform(int f[])
73 throws IllegalArgumentException {
74 return fht(f);
75 }
76
77 /**
78 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
79 * <br>
80 * Requires <b>Nlog2N = n2</b><sup>n</sup> additions.
81 * <br>
82 * <br>
83 * <b><u>Short Table of manual calculation for N=8:</u></b>
84 * <ol>
85 * <li><b>x</b> is the input vector we want to transform</li>
86 * <li><b>y</b> is the output vector which is our desired result</li>
87 * <li>a and b are just helper rows</li>
88 * </ol>
89 * <pre>
90 * <code>
91 * +----+----------+---------+----------+
92 * | <b>x</b> | <b>a</b> | <b>b</b> | <b>y</b> |
93 * +----+----------+---------+----------+
94 * | x<sub>0</sub> | a<sub>0</sub>=x<sub>0</sub>+x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>+a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>+b<sub>1</sub> |
95 * +----+----------+---------+----------+
96 * | x<sub>1</sub> | a<sub>1</sub>=x<sub>2</sub>+x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>+a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>+b<sub>3</sub> |
97 * +----+----------+---------+----------+
98 * | x<sub>2</sub> | a<sub>2</sub>=x<sub>4</sub>+x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>+a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>+b<sub>5</sub> |
99 * +----+----------+---------+----------+
100 * | x<sub>3</sub> | a<sub>3</sub>=x<sub>6</sub>+x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>+a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>+b<sub>7</sub> |
101 * +----+----------+---------+----------+
102 * | x<sub>4</sub> | a<sub>0</sub>=x<sub>0</sub>-x<sub>1</sub> | b<sub>0</sub>=a<sub>0</sub>-a<sub>1</sub> | y<sub>0</sub>=b<sub>0</sub>-b<sub>1</sub> |
103 * +----+----------+---------+----------+
104 * | x<sub>5</sub> | a<sub>1</sub>=x<sub>2</sub>-x<sub>3</sub> | b<sub>0</sub>=a<sub>2</sub>-a<sub>3</sub> | y<sub>0</sub>=b<sub>2</sub>-b<sub>3</sub> |
105 * +----+----------+---------+----------+
106 * | x<sub>6</sub> | a<sub>2</sub>=x<sub>4</sub>-x<sub>5</sub> | b<sub>0</sub>=a<sub>4</sub>-a<sub>5</sub> | y<sub>0</sub>=b<sub>4</sub>-b<sub>5</sub> |
107 * +----+----------+---------+----------+
108 * | x<sub>7</sub> | a<sub>3</sub>=x<sub>6</sub>-x<sub>7</sub> | b<sub>0</sub>=a<sub>6</sub>-a<sub>7</sub> | y<sub>0</sub>=b<sub>6</sub>-b<sub>7</sub> |
109 * +----+----------+---------+----------+
110 * </code>
111 * </pre>
112 *
113 * <b><u>How it works</u></b>
114 * <ol>
115 * <li>Construct a matrix with N rows and n+1 columns<br> <b>hadm[n+1][N]</b>
116 * <br><i>(If I use [x][y] it always means [row-offset][column-offset] of a Matrix with n rows and m columns. Its entries go from M[0][0] to M[n][m])</i></li>
117 * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li>
118 * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows.
119 * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1].
120 * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column
121 * </li>
122 * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows.
123 * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1].
124 * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column
125 * </li>
126 * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above.
127 * <li>The output vector y is now in the last column of <b>hadm</b></li>
128 * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li>
129 * </ol>
130 * <br>
131 * <b><u>Visually</u></b>
132 * <pre>
133 * +--------+---+---+---+-----+---+
134 * | 0 | 1 | 2 | 3 | ... |n+1|
135 * +------+--------+---+---+---+-----+---+
136 * |0 | x<sub>0</sub> | /\ |
137 * |1 | x<sub>1</sub> | || |
138 * |2 | x<sub>2</sub> | <= D<sub>top</sub> => |
139 * |... | ... | || |
140 * |N/2-1 | x<sub>N/2-1</sub> | \/ |
141 * +------+--------+---+---+---+-----+---+
142 * |N/2 | x<sub>N/2</sub> | /\ |
143 * |N/2+1 | x<sub>N/2+1</sub> | || |
144 * |N/2+2 | x<sub>N/2+2</sub> | <= D<sub>bottom</sub> => | which is in the last column of the matrix
145 * |... | ... | || |
146 * |N | x<sub>N/2</sub> | \/ |
147 * +------+--------+---+---+---+-----+---+
148 * </pre>
149 *
150 * @param x input vector
151 * @return y output vector
152 * @exception IllegalArgumentException if input array is not a power of 2
153 */
154 protected double[] fht(double x[]) throws IllegalArgumentException {
155
156 // n is the row count of the input vector x
157 final int n = x.length;
158 final int halfN = n / 2;
159
160 // n has to be of the form n = 2^p !!
161 if (!FastFourierTransformer.isPowerOf2(n)) {
162 throw MathRuntimeException.createIllegalArgumentException(
163 "{0} is not a power of 2",
164 n);
165 }
166
167 // Instead of creating a matrix with p+1 columns and n rows
168 // we will use two single dimension arrays which we will use in an alternating way.
169 double[] yPrevious = new double[n];
170 double[] yCurrent = x.clone();
171
172 // iterate from left to right (column)
173 for (int j = 1; j < n; j <<= 1) {
174
175 // switch columns
176 final double[] yTmp = yCurrent;
177 yCurrent = yPrevious;
178 yPrevious = yTmp;
179
180 // iterate from top to bottom (row)
181 for (int i = 0; i < halfN; ++i) {
182 // D<sub>top</sub>
183 // The top part works with addition
184 final int twoI = 2 * i;
185 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
186 }
187 for (int i = halfN; i < n; ++i) {
188 // D<sub>bottom</sub>
189 // The bottom part works with subtraction
190 final int twoI = 2 * i;
191 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
192 }
193 }
194
195 // return the last computed output vector y
196 return yCurrent;
197
198 }
199 /**
200 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
201 * @param x input vector
202 * @return y output vector
203 * @exception IllegalArgumentException if input array is not a power of 2
204 */
205 protected int[] fht(int x[]) throws IllegalArgumentException {
206
207 // n is the row count of the input vector x
208 final int n = x.length;
209 final int halfN = n / 2;
210
211 // n has to be of the form n = 2^p !!
212 if (!FastFourierTransformer.isPowerOf2(n)) {
213 throw MathRuntimeException.createIllegalArgumentException(
214 "{0} is not a power of 2",
215 n);
216 }
217
218 // Instead of creating a matrix with p+1 columns and n rows
219 // we will use two single dimension arrays which we will use in an alternating way.
220 int[] yPrevious = new int[n];
221 int[] yCurrent = x.clone();
222
223 // iterate from left to right (column)
224 for (int j = 1; j < n; j <<= 1) {
225
226 // switch columns
227 final int[] yTmp = yCurrent;
228 yCurrent = yPrevious;
229 yPrevious = yTmp;
230
231 // iterate from top to bottom (row)
232 for (int i = 0; i < halfN; ++i) {
233 // D<sub>top</sub>
234 // The top part works with addition
235 final int twoI = 2 * i;
236 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
237 }
238 for (int i = halfN; i < n; ++i) {
239 // D<sub>bottom</sub>
240 // The bottom part works with subtraction
241 final int twoI = 2 * i;
242 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
243 }
244 }
245
246 // return the last computed output vector y
247 return yCurrent;
248
249 }
250
251 }