View Javadoc

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.analysis.*;
20  import org.apache.commons.math.complex.*;
21  import org.apache.commons.math.FunctionEvaluationException;
22  import org.apache.commons.math.MathRuntimeException;
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   * <p>As of version 2.0 this no longer implements Serializable</p>
37   * 
38   * @version $Revision: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 2009) $
39   * @since 1.2
40   */
41  public class FastSineTransformer implements RealTransformer {
42  
43      /**
44       * Construct a default transformer.
45       */
46      public FastSineTransformer() {
47          super();
48      }
49  
50      /**
51       * Transform the given real data set.
52       * <p>
53       * The formula is F<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
54       * </p>
55       * 
56       * @param f the real data array to be transformed
57       * @return the real transformed array
58       * @throws IllegalArgumentException if any parameters are invalid
59       */
60      public double[] transform(double f[])
61          throws IllegalArgumentException {
62          return fst(f);
63      }
64  
65      /**
66       * Transform the given real function, sampled on the given interval.
67       * <p>
68       * The formula is F<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
69       * </p>
70       * 
71       * @param f the function to be sampled and transformed
72       * @param min the lower bound for the interval
73       * @param max the upper bound for the interval
74       * @param n the number of sample points
75       * @return the real transformed array
76       * @throws FunctionEvaluationException if function cannot be evaluated
77       * at some point
78       * @throws IllegalArgumentException if any parameters are invalid
79       */
80      public double[] transform(UnivariateRealFunction f,
81                                double min, double max, int n)
82          throws FunctionEvaluationException, IllegalArgumentException {
83  
84          double data[] = FastFourierTransformer.sample(f, min, max, n);
85          data[0] = 0.0;
86          return fst(data);
87      }
88  
89      /**
90       * Transform the given real data set.
91       * <p>
92       * The formula is F<sub>n</sub> = &radic;(2/N) &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
93       * </p>
94       * 
95       * @param f the real data array to be transformed
96       * @return the real transformed array
97       * @throws IllegalArgumentException if any parameters are invalid
98       */
99      public double[] transform2(double f[]) throws IllegalArgumentException {
100 
101         double scaling_coefficient = Math.sqrt(2.0 / f.length);
102         return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
103     }
104 
105     /**
106      * Transform the given real function, sampled on the given interval.
107      * <p>
108      * The formula is F<sub>n</sub> = &radic;(2/N) &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
109      * </p>
110      * 
111      * @param f the function to be sampled and transformed
112      * @param min the lower bound for the interval
113      * @param max the upper bound for the interval
114      * @param n the number of sample points
115      * @return the real transformed array
116      * @throws FunctionEvaluationException if function cannot be evaluated
117      * at some point
118      * @throws IllegalArgumentException if any parameters are invalid
119      */
120     public double[] transform2(
121         UnivariateRealFunction f, double min, double max, int n)
122         throws FunctionEvaluationException, IllegalArgumentException {
123 
124         double data[] = FastFourierTransformer.sample(f, min, max, n);
125         data[0] = 0.0;
126         double scaling_coefficient = Math.sqrt(2.0 / n);
127         return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
128     }
129 
130     /**
131      * Inversely transform the given real data set.
132      * <p>
133      * The formula is f<sub>k</sub> = (2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
134      * </p>
135      * 
136      * @param f the real data array to be inversely transformed
137      * @return the real inversely transformed array
138      * @throws IllegalArgumentException if any parameters are invalid
139      */
140     public double[] inversetransform(double f[]) throws IllegalArgumentException {
141 
142         double scaling_coefficient = 2.0 / f.length;
143         return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
144     }
145 
146     /**
147      * Inversely transform the given real function, sampled on the given interval.
148      * <p>
149      * The formula is f<sub>k</sub> = (2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
150      * </p>
151      * 
152      * @param f the function to be sampled and inversely transformed
153      * @param min the lower bound for the interval
154      * @param max the upper bound for the interval
155      * @param n the number of sample points
156      * @return the real inversely transformed array
157      * @throws FunctionEvaluationException if function cannot be evaluated
158      * at some point
159      * @throws IllegalArgumentException if any parameters are invalid
160      */
161     public double[] inversetransform(UnivariateRealFunction f, double min, double max, int n)
162         throws FunctionEvaluationException, IllegalArgumentException {
163 
164         double data[] = FastFourierTransformer.sample(f, min, max, n);
165         data[0] = 0.0;
166         double scaling_coefficient = 2.0 / n;
167         return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
168     }
169 
170     /**
171      * Inversely transform the given real data set.
172      * <p>
173      * The formula is f<sub>k</sub> = &radic;(2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
174      * </p>
175      * 
176      * @param f the real data array to be inversely transformed
177      * @return the real inversely transformed array
178      * @throws IllegalArgumentException if any parameters are invalid
179      */
180     public double[] inversetransform2(double f[]) throws IllegalArgumentException {
181 
182         return transform2(f);
183     }
184 
185     /**
186      * Inversely transform the given real function, sampled on the given interval.
187      * <p>
188      * The formula is f<sub>k</sub> = &radic;(2/N) &sum;<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(&pi; nk/N)
189      * </p>
190      * 
191      * @param f the function to be sampled and inversely transformed
192      * @param min the lower bound for the interval
193      * @param max the upper bound for the interval
194      * @param n the number of sample points
195      * @return the real inversely transformed array
196      * @throws FunctionEvaluationException if function cannot be evaluated
197      * at some point
198      * @throws IllegalArgumentException if any parameters are invalid
199      */
200     public double[] inversetransform2(UnivariateRealFunction f, double min, double max, int n)
201         throws FunctionEvaluationException, IllegalArgumentException {
202 
203         return transform2(f, min, max, n);
204     }
205 
206     /**
207      * Perform the FST algorithm (including inverse).
208      *
209      * @param f the real data array to be transformed
210      * @return the real transformed array
211      * @throws IllegalArgumentException if any parameters are invalid
212      */
213     protected double[] fst(double f[]) throws IllegalArgumentException {
214 
215         double A, B, x[], F[] = new double[f.length];
216 
217         FastFourierTransformer.verifyDataSet(f);
218         if (f[0] != 0.0) {
219             throw MathRuntimeException.createIllegalArgumentException(
220                     "first element is not 0: {0}",
221                     f[0]);
222         }
223         int N = f.length;
224         if (N == 1) {       // trivial case
225             F[0] = 0.0;
226             return F;
227         }
228 
229         // construct a new array and perform FFT on it
230         x = new double[N];
231         x[0] = 0.0;
232         x[N >> 1] = 2.0 * f[N >> 1];
233         for (int i = 1; i < (N >> 1); i++) {
234             A = Math.sin(i * Math.PI / N) * (f[i] + f[N-i]);
235             B = 0.5 * (f[i] - f[N-i]);
236             x[i] = A + B;
237             x[N-i] = A - B;
238         }
239         FastFourierTransformer transformer = new FastFourierTransformer();
240         Complex y[] = transformer.transform(x);
241 
242         // reconstruct the FST result for the original array
243         F[0] = 0.0;
244         F[1] = 0.5 * y[0].getReal();
245         for (int i = 1; i < (N >> 1); i++) {
246             F[2*i] = -y[i].getImaginary();
247             F[2*i+1] = y[i].getReal() + F[2*i-1];
248         }
249 
250         return F;
251     }
252 }