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 */
17package org.apache.commons.math.transform;
18
19import org.apache.commons.math.FunctionEvaluationException;
20import org.apache.commons.math.MathRuntimeException;
21import org.apache.commons.math.analysis.UnivariateRealFunction;
22import org.apache.commons.math.exception.util.LocalizedFormats;
23
24/**
25 * Implements the <a href="http://www.archive.chipcenter.com/dsp/DSP000517F1.html">Fast Hadamard Transform</a> (FHT).
26 * Transformation of an input vector x to the output vector y.
27 * <p>In addition to transformation of real vectors, the Hadamard transform can
28 * transform integer vectors into integer vectors. However, this integer transform
29 * cannot be inverted directly. Due to a scaling factor it may lead to rational results.
30 * As an example, the inverse transform of integer vector (0, 1, 0, 1) is rational
31 * vector (1/2, -1/2, 0, 0).</p>
32 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
33 * @since 2.0
34 */
35public class FastHadamardTransformer implements RealTransformer {
36
37    /** {@inheritDoc} */
38    public double[] transform(double f[])
39        throws IllegalArgumentException {
40        return fht(f);
41    }
42
43    /** {@inheritDoc} */
44    public double[] transform(UnivariateRealFunction f,
45                              double min, double max, int n)
46        throws FunctionEvaluationException, IllegalArgumentException {
47        return fht(FastFourierTransformer.sample(f, min, max, n));
48    }
49
50    /** {@inheritDoc} */
51    public double[] inversetransform(double f[])
52    throws IllegalArgumentException {
53        return FastFourierTransformer.scaleArray(fht(f), 1.0 / f.length);
54   }
55
56    /** {@inheritDoc} */
57    public double[] inversetransform(UnivariateRealFunction f,
58                                     double min, double max, int n)
59        throws FunctionEvaluationException, IllegalArgumentException {
60        final double[] unscaled =
61            fht(FastFourierTransformer.sample(f, min, max, n));
62        return FastFourierTransformer.scaleArray(unscaled, 1.0 / n);
63    }
64
65    /**
66     * Transform the given real data set.
67     * <p>The integer transform cannot be inverted directly, due to a scaling
68     * factor it may lead to double results.</p>
69     * @param f the integer data array to be transformed (signal)
70     * @return the integer transformed array (spectrum)
71     * @throws IllegalArgumentException if any parameters are invalid
72     */
73    public int[] transform(int f[])
74        throws IllegalArgumentException {
75        return fht(f);
76    }
77
78    /**
79     * The FHT (Fast Hadamard Transformation) which uses only subtraction and addition.
80     * <br>
81     * Requires <b>Nlog2N = n2</b><sup>n</sup> additions.
82     * <br>
83     * <br>
84     * <b><u>Short Table of manual calculation for N=8:</u></b>
85     * <ol>
86     * <li><b>x</b> is the input vector we want to transform</li>
87     * <li><b>y</b> is the output vector which is our desired result</li>
88     * <li>a and b are just helper rows</li>
89     * </ol>
90     * <pre>
91     * <code>
92     * +----+----------+---------+----------+
93     * | <b>x</b>  |    <b>a</b>     |    <b>b</b>    |    <b>y</b>     |
94     * +----+----------+---------+----------+
95     * | 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> |
96     * +----+----------+---------+----------+
97     * | 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> |
98     * +----+----------+---------+----------+
99     * | 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}
253