~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/java/org/apache/commons/math/transform/FastSineTransformer.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
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 java.io.Serializable;
20
 
import org.apache.commons.math.analysis.*;
21
 
import org.apache.commons.math.complex.*;
22
 
import org.apache.commons.math.MathException;
23
 
 
24
 
/**
25
 
 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
26
 
 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a>
27
 
 * for transformation of one-dimensional data sets. For reference, see
28
 
 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
29
 
 * <p>
30
 
 * FST is its own inverse, up to a multiplier depending on conventions.
31
 
 * The equations are listed in the comments of the corresponding methods.</p>
32
 
 * <p>
33
 
 * Similar to FFT, we also require the length of data set to be power of 2.
34
 
 * In addition, the first element must be 0 and it's enforced in function
35
 
 * transformation after sampling.</p>
36
 
 *
37
 
 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
38
 
 * @since 1.2
39
 
 */
40
 
public class FastSineTransformer implements Serializable {
41
 
 
42
 
    /** serializable version identifier */
43
 
    static final long serialVersionUID = -478002039949390854L;
44
 
 
45
 
    /**
46
 
     * Construct a default transformer.
47
 
     */
48
 
    public FastSineTransformer() {
49
 
        super();
50
 
    }
51
 
 
52
 
    /**
53
 
     * Transform the given real data set.
54
 
     * <p>
55
 
     * The formula is $ F_n = \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
56
 
     * </p>
57
 
     * 
58
 
     * @param f the real data array to be transformed
59
 
     * @return the real transformed array
60
 
     * @throws MathException if any math-related errors occur
61
 
     * @throws IllegalArgumentException if any parameters are invalid
62
 
     */
63
 
    public double[] transform(double f[]) throws MathException,
64
 
        IllegalArgumentException {
65
 
 
66
 
        return fst(f);
67
 
    }
68
 
 
69
 
    /**
70
 
     * Transform the given real function, sampled on the given interval.
71
 
     * <p>
72
 
     * The formula is $ F_n = \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
73
 
     * </p>
74
 
     * 
75
 
     * @param f the function to be sampled and transformed
76
 
     * @param min the lower bound for the interval
77
 
     * @param max the upper bound for the interval
78
 
     * @param n the number of sample points
79
 
     * @return the real transformed array
80
 
     * @throws MathException if any math-related errors occur
81
 
     * @throws IllegalArgumentException if any parameters are invalid
82
 
     */
83
 
    public double[] transform(
84
 
        UnivariateRealFunction f, double min, double max, int n)
85
 
        throws MathException, IllegalArgumentException {
86
 
 
87
 
        double data[] = FastFourierTransformer.sample(f, min, max, n);
88
 
        data[0] = 0.0;
89
 
        return fst(data);
90
 
    }
91
 
 
92
 
    /**
93
 
     * Transform the given real data set.
94
 
     * <p>
95
 
     * The formula is $ F_n = \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
96
 
     * </p>
97
 
     * 
98
 
     * @param f the real data array to be transformed
99
 
     * @return the real transformed array
100
 
     * @throws MathException if any math-related errors occur
101
 
     * @throws IllegalArgumentException if any parameters are invalid
102
 
     */
103
 
    public double[] transform2(double f[]) throws MathException,
104
 
        IllegalArgumentException {
105
 
 
106
 
        double scaling_coefficient = Math.sqrt(2.0 / f.length);
107
 
        return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
108
 
    }
109
 
 
110
 
    /**
111
 
     * Transform the given real function, sampled on the given interval.
112
 
     * <p>
113
 
     * The formula is $ F_n = \sqrt{2/N} \Sigma_{k=0}^{N-1} f_k \sin(\pi nk/N) $
114
 
     * </p>
115
 
     * 
116
 
     * @param f the function to be sampled and transformed
117
 
     * @param min the lower bound for the interval
118
 
     * @param max the upper bound for the interval
119
 
     * @param n the number of sample points
120
 
     * @return the real transformed array
121
 
     * @throws MathException if any math-related errors occur
122
 
     * @throws IllegalArgumentException if any parameters are invalid
123
 
     */
124
 
    public double[] transform2(
125
 
        UnivariateRealFunction f, double min, double max, int n)
126
 
        throws MathException, IllegalArgumentException {
127
 
 
128
 
        double data[] = FastFourierTransformer.sample(f, min, max, n);
129
 
        data[0] = 0.0;
130
 
        double scaling_coefficient = Math.sqrt(2.0 / n);
131
 
        return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
132
 
    }
133
 
 
134
 
    /**
135
 
     * Inversely transform the given real data set.
136
 
     * <p>
137
 
     * The formula is $ f_k = (2/N) \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
138
 
     * </p>
139
 
     * 
140
 
     * @param f the real data array to be inversely transformed
141
 
     * @return the real inversely transformed array
142
 
     * @throws MathException if any math-related errors occur
143
 
     * @throws IllegalArgumentException if any parameters are invalid
144
 
     */
145
 
    public double[] inversetransform(double f[]) throws MathException,
146
 
        IllegalArgumentException {
147
 
 
148
 
        double scaling_coefficient = 2.0 / f.length;
149
 
        return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
150
 
    }
151
 
 
152
 
    /**
153
 
     * Inversely transform the given real function, sampled on the given interval.
154
 
     * <p>
155
 
     * The formula is $ f_k = (2/N) \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
156
 
     * </p>
157
 
     * 
158
 
     * @param f the function to be sampled and inversely transformed
159
 
     * @param min the lower bound for the interval
160
 
     * @param max the upper bound for the interval
161
 
     * @param n the number of sample points
162
 
     * @return the real inversely transformed array
163
 
     * @throws MathException if any math-related errors occur
164
 
     * @throws IllegalArgumentException if any parameters are invalid
165
 
     */
166
 
    public double[] inversetransform(
167
 
        UnivariateRealFunction f, double min, double max, int n)
168
 
        throws MathException, IllegalArgumentException {
169
 
 
170
 
        double data[] = FastFourierTransformer.sample(f, min, max, n);
171
 
        data[0] = 0.0;
172
 
        double scaling_coefficient = 2.0 / n;
173
 
        return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
174
 
    }
175
 
 
176
 
    /**
177
 
     * Inversely transform the given real data set.
178
 
     * <p>
179
 
     * The formula is $ f_k = \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
180
 
     * </p>
181
 
     * 
182
 
     * @param f the real data array to be inversely transformed
183
 
     * @return the real inversely transformed array
184
 
     * @throws MathException if any math-related errors occur
185
 
     * @throws IllegalArgumentException if any parameters are invalid
186
 
     */
187
 
    public double[] inversetransform2(double f[]) throws MathException,
188
 
        IllegalArgumentException {
189
 
 
190
 
        return transform2(f);
191
 
    }
192
 
 
193
 
    /**
194
 
     * Inversely transform the given real function, sampled on the given interval.
195
 
     * <p>
196
 
     * The formula is $ f_k = \sqrt{2/N} \Sigma_{n=0}^{N-1} F_n \sin(\pi nk/N) $
197
 
     * </p>
198
 
     * 
199
 
     * @param f the function to be sampled and inversely transformed
200
 
     * @param min the lower bound for the interval
201
 
     * @param max the upper bound for the interval
202
 
     * @param n the number of sample points
203
 
     * @return the real inversely transformed array
204
 
     * @throws MathException if any math-related errors occur
205
 
     * @throws IllegalArgumentException if any parameters are invalid
206
 
     */
207
 
    public double[] inversetransform2(
208
 
        UnivariateRealFunction f, double min, double max, int n)
209
 
        throws MathException, IllegalArgumentException {
210
 
 
211
 
        return transform2(f, min, max, n);
212
 
    }
213
 
 
214
 
    /**
215
 
     * Perform the FST algorithm (including inverse).
216
 
     *
217
 
     * @param f the real data array to be transformed
218
 
     * @return the real transformed array
219
 
     * @throws MathException if any math-related errors occur
220
 
     * @throws IllegalArgumentException if any parameters are invalid
221
 
     */
222
 
    protected double[] fst(double f[]) throws MathException,
223
 
        IllegalArgumentException {
224
 
 
225
 
        double A, B, x[], F[] = new double[f.length];
226
 
 
227
 
        FastFourierTransformer.verifyDataSet(f);
228
 
        if (f[0] != 0.0) {
229
 
            throw new IllegalArgumentException
230
 
                ("The first element is not zero: " + f[0]);
231
 
        }
232
 
        int N = f.length;
233
 
        if (N == 1) {       // trivial case
234
 
            F[0] = 0.0;
235
 
            return F;
236
 
        }
237
 
 
238
 
        // construct a new array and perform FFT on it
239
 
        x = new double[N];
240
 
        x[0] = 0.0;
241
 
        x[N >> 1] = 2.0 * f[N >> 1];
242
 
        for (int i = 1; i < (N >> 1); i++) {
243
 
            A = Math.sin(i * Math.PI / N) * (f[i] + f[N-i]);
244
 
            B = 0.5 * (f[i] - f[N-i]);
245
 
            x[i] = A + B;
246
 
            x[N-i] = A - B;
247
 
        }
248
 
        FastFourierTransformer transformer = new FastFourierTransformer();
249
 
        Complex y[] = transformer.transform(x);
250
 
 
251
 
        // reconstruct the FST result for the original array
252
 
        F[0] = 0.0;
253
 
        F[1] = 0.5 * y[0].getReal();
254
 
        for (int i = 1; i < (N >> 1); i++) {
255
 
            F[2*i] = -y[i].getImaginary();
256
 
            F[2*i+1] = y[i].getReal() + F[2*i-1];
257
 
        }
258
 
 
259
 
        return F;
260
 
    }
261
 
}