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.math.transform;
018
019 import org.apache.commons.math.FunctionEvaluationException;
020 import org.apache.commons.math.MathRuntimeException;
021 import org.apache.commons.math.analysis.UnivariateRealFunction;
022 import org.apache.commons.math.exception.util.LocalizedFormats;
023
024 /**
025 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
026 * Transformation of an input vector x to the output vector y.
027 * <p>In addition to transformation of real vectors, the Hadamard transform can
028 * transform integer vectors into integer vectors. However, this integer transform
029 * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
030 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
031 * vector (1/2, -1/2, 0, 0).</p>
032 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
033 * @since 2.0
034 */
035 public class FastHadamardTransformer implements RealTransformer {
036
037 /** {@inheritDoc} */
038 public double[] transform(double f[])
039 throws IllegalArgumentException {
040 return fht(f);
041 }
042
043 /** {@inheritDoc} */
044 public double[] transform(UnivariateRealFunction f,
045 double min, double max, int n)
046 throws FunctionEvaluationException, IllegalArgumentException {
047 return fht(FastFourierTransformer.sample(f, min, max, n));
048 }
049
050 /** {@inheritDoc} */
051 public double[] inversetransform(double f[])
052 throws IllegalArgumentException {
053 return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
054 }
055
056 /** {@inheritDoc} */
057 public double[] inversetransform(UnivariateRealFunction f,
058 double min, double max, int n)
059 throws FunctionEvaluationException, IllegalArgumentException {
060 final double[] unscaled =
061 fht(FastFourierTransformer.sample(f, min, max, n));
062 return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
063 }
064
065 /**
066 * Transform the given real data set.
067 * <p>The integer transform cannot be inverted directly, due to a scaling
068 * factor it may lead to double results.</p>
069 * @param f the integer data array to be transformed (signal)
070 * @return the integer transformed array (spectrum)
071 * @throws IllegalArgumentException if any parameters are invalid
072 */
073 public int[] transform(int f[])
074 throws IllegalArgumentException {
075 return fht(f);
076 }
077
078 /**
079 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
080 * <br>
081 * Requires <b>Nlog2N = n2</b><sup>n</sup> additions.
082 * <br>
083 * <br>
084 * <b><u>Short Table of manual calculation for N=8:</u></b>
085 * <ol>
086 * <li><b>x</b> is the input vector we want to transform</li>
087 * <li><b>y</b> is the output vector which is our desired result</li>
088 * <li>a and b are just helper rows</li>
089 * </ol>
090 * <pre>
091 * <code>
092 * +----+----------+---------+----------+
093 * | <b>x</b> | <b>a</b> | <b>b</b> | <b>y</b> |
094 * +----+----------+---------+----------+
095 * | 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> |
096 * +----+----------+---------+----------+
097 * | 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> |
098 * +----+----------+---------+----------+
099 * | 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> |
100 * +----+----------+---------+----------+
101 * | 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> |
102 * +----+----------+---------+----------+
103 * | 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> |
104 * +----+----------+---------+----------+
105 * | 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> |
106 * +----+----------+---------+----------+
107 * | 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> |
108 * +----+----------+---------+----------+
109 * | 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> |
110 * +----+----------+---------+----------+
111 * </code>
112 * </pre>
113 *
114 * <b><u>How it works</u></b>
115 * <ol>
116 * <li>Construct a matrix with N rows and n+1 columns<br> <b>hadm[n+1][N]</b>
117 * <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>
118 * <li>Place the input vector <b>x[N]</b> in the first column of the matrix <b>hadm</b></li>
119 * <li>The entries of the submatrix D<sub>top</sub> are calculated as follows.
120 * <br>D<sub>top</sub> goes from entry [0][1] to [N/2-1][n+1].
121 * <br>The columns of D<sub>top</sub> are the pairwise mutually exclusive sums of the previous column
122 * </li>
123 * <li>The entries of the submatrix D<sub>bottom</sub> are calculated as follows.
124 * <br>D<sub>bottom</sub> goes from entry [N/2][1] to [N][n+1].
125 * <br>The columns of D<sub>bottom</sub> are the pairwise differences of the previous column
126 * </li>
127 * <li>How D<sub>top</sub> and D<sub>bottom</sub> you can understand best with the example for N=8 above.
128 * <li>The output vector y is now in the last column of <b>hadm</b></li>
129 * <li><i>Algorithm from: http://www.archive.chipcenter.com/dsp/DSP000517F1.html</i></li>
130 * </ol>
131 * <br>
132 * <b><u>Visually</u></b>
133 * <pre>
134 * +--------+---+---+---+-----+---+
135 * | 0 | 1 | 2 | 3 | ... |n+1|
136 * +------+--------+---+---+---+-----+---+
137 * |0 | x<sub>0</sub> | /\ |
138 * |1 | x<sub>1</sub> | || |
139 * |2 | x<sub>2</sub> | <= D<sub>top</sub> => |
140 * |... | ... | || |
141 * |N/2-1 | x<sub>N/2-1</sub> | \/ |
142 * +------+--------+---+---+---+-----+---+
143 * |N/2 | x<sub>N/2</sub> | /\ |
144 * |N/2+1 | x<sub>N/2+1</sub> | || |
145 * |N/2+2 | x<sub>N/2+2</sub> | <= D<sub>bottom</sub> => | which is in the last column of the matrix
146 * |... | ... | || |
147 * |N | x<sub>N/2</sub> | \/ |
148 * +------+--------+---+---+---+-----+---+
149 * </pre>
150 *
151 * @param x input vector
152 * @return y output vector
153 * @exception IllegalArgumentException if input array is not a power of 2
154 */
155 protected double[] fht(double x[]) throws IllegalArgumentException {
156
157 // n is the row count of the input vector x
158 final int n = x.length;
159 final int halfN = n / 2;
160
161 // n has to be of the form n = 2^p !!
162 if (!FastFourierTransformer.isPowerOf2(n)) {
163 throw MathRuntimeException.createIllegalArgumentException(
164 LocalizedFormats.NOT_POWER_OF_TWO,
165 n);
166 }
167
168 // Instead of creating a matrix with p+1 columns and n rows
169 // we will use two single dimension arrays which we will use in an alternating way.
170 double[] yPrevious = new double[n];
171 double[] yCurrent = x.clone();
172
173 // iterate from left to right (column)
174 for (int j = 1; j < n; j <<= 1) {
175
176 // switch columns
177 final double[] yTmp = yCurrent;
178 yCurrent = yPrevious;
179 yPrevious = yTmp;
180
181 // iterate from top to bottom (row)
182 for (int i = 0; i < halfN; ++i) {
183 // D<sub>top</sub>
184 // The top part works with addition
185 final int twoI = 2 * i;
186 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
187 }
188 for (int i = halfN; i < n; ++i) {
189 // D<sub>bottom</sub>
190 // The bottom part works with subtraction
191 final int twoI = 2 * i;
192 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
193 }
194 }
195
196 // return the last computed output vector y
197 return yCurrent;
198
199 }
200 /**
201 * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
202 * @param x input vector
203 * @return y output vector
204 * @exception IllegalArgumentException if input array is not a power of 2
205 */
206 protected int[] fht(int x[]) throws IllegalArgumentException {
207
208 // n is the row count of the input vector x
209 final int n = x.length;
210 final int halfN = n / 2;
211
212 // n has to be of the form n = 2^p !!
213 if (!FastFourierTransformer.isPowerOf2(n)) {
214 throw MathRuntimeException.createIllegalArgumentException(
215 LocalizedFormats.NOT_POWER_OF_TWO,
216 n);
217 }
218
219 // Instead of creating a matrix with p+1 columns and n rows
220 // we will use two single dimension arrays which we will use in an alternating way.
221 int[] yPrevious = new int[n];
222 int[] yCurrent = x.clone();
223
224 // iterate from left to right (column)
225 for (int j = 1; j < n; j <<= 1) {
226
227 // switch columns
228 final int[] yTmp = yCurrent;
229 yCurrent = yPrevious;
230 yPrevious = yTmp;
231
232 // iterate from top to bottom (row)
233 for (int i = 0; i < halfN; ++i) {
234 // D<sub>top</sub>
235 // The top part works with addition
236 final int twoI = 2 * i;
237 yCurrent[i] = yPrevious[twoI] + yPrevious[twoI + 1];
238 }
239 for (int i = halfN; i < n; ++i) {
240 // D<sub>bottom</sub>
241 // The bottom part works with subtraction
242 final int twoI = 2 * i;
243 yCurrent[i] = yPrevious[twoI - n] - yPrevious[twoI - n + 1];
244 }
245 }
246
247 // return the last computed output vector y
248 return yCurrent;
249
250 }
251
252 }