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.complex.Complex;
023
024 /**
025 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
026 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a>
027 * for transformation of one-dimensional data sets. For reference, see
028 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
029 * <p>
030 * FST is its own inverse, up to a multiplier depending on conventions.
031 * The equations are listed in the comments of the corresponding methods.</p>
032 * <p>
033 * Similar to FFT, we also require the length of data set to be power of 2.
034 * In addition, the first element must be 0 and it's enforced in function
035 * transformation after sampling.</p>
036 * <p>As of version 2.0 this no longer implements Serializable</p>
037 *
038 * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
039 * @since 1.2
040 */
041 public class FastSineTransformer implements RealTransformer {
042
043 /**
044 * Construct a default transformer.
045 */
046 public FastSineTransformer() {
047 super();
048 }
049
050 /**
051 * Transform the given real data set.
052 * <p>
053 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N)
054 * </p>
055 *
056 * @param f the real data array to be transformed
057 * @return the real transformed array
058 * @throws IllegalArgumentException if any parameters are invalid
059 */
060 public double[] transform(double f[])
061 throws IllegalArgumentException {
062 return fst(f);
063 }
064
065 /**
066 * Transform the given real function, sampled on the given interval.
067 * <p>
068 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N)
069 * </p>
070 *
071 * @param f the function to be sampled and transformed
072 * @param min the lower bound for the interval
073 * @param max the upper bound for the interval
074 * @param n the number of sample points
075 * @return the real transformed array
076 * @throws FunctionEvaluationException if function cannot be evaluated
077 * at some point
078 * @throws IllegalArgumentException if any parameters are invalid
079 */
080 public double[] transform(UnivariateRealFunction f,
081 double min, double max, int n)
082 throws FunctionEvaluationException, IllegalArgumentException {
083
084 double data[] = FastFourierTransformer.sample(f, min, max, n);
085 data[0] = 0.0;
086 return fst(data);
087 }
088
089 /**
090 * Transform the given real data set.
091 * <p>
092 * The formula is F<sub>n</sub> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N)
093 * </p>
094 *
095 * @param f the real data array to be transformed
096 * @return the real transformed array
097 * @throws IllegalArgumentException if any parameters are invalid
098 */
099 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 final double transformed[] = 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 final int n = f.length;
224 if (n == 1) { // trivial case
225 transformed[0] = 0.0;
226 return transformed;
227 }
228
229 // construct a new array and perform FFT on it
230 final double[] 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 final double a = Math.sin(i * Math.PI / n) * (f[i] + f[n-i]);
235 final double 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 transformed[0] = 0.0;
244 transformed[1] = 0.5 * y[0].getReal();
245 for (int i = 1; i < (n >> 1); i++) {
246 transformed[2 * i] = -y[i].getImaginary();
247 transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1];
248 }
249
250 return transformed;
251 }
252 }