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> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π 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> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π 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> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π 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> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π 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) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π 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) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π 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> = √(2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π 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> = √(2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π 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 }