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    import org.apache.commons.math.exception.util.LocalizedFormats;
024    import org.apache.commons.math.util.FastMath;
025    
026    /**
027     * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
028     * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a>
029     * for transformation of one-dimensional data sets. For reference, see
030     * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3.
031     * <p>
032     * FST is its own inverse, up to a multiplier depending on conventions.
033     * The equations are listed in the comments of the corresponding methods.</p>
034     * <p>
035     * Similar to FFT, we also require the length of data set to be power of 2.
036     * In addition, the first element must be 0 and it's enforced in function
037     * transformation after sampling.</p>
038     * <p>As of version 2.0 this no longer implements Serializable</p>
039     *
040     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
041     * @since 1.2
042     */
043    public class FastSineTransformer implements RealTransformer {
044    
045        /**
046         * Construct a default transformer.
047         */
048        public FastSineTransformer() {
049            super();
050        }
051    
052        /**
053         * Transform the given real data set.
054         * <p>
055         * The formula is F<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
056         * </p>
057         *
058         * @param f the real data array to be transformed
059         * @return the real transformed array
060         * @throws IllegalArgumentException if any parameters are invalid
061         */
062        public double[] transform(double f[])
063            throws IllegalArgumentException {
064            return fst(f);
065        }
066    
067        /**
068         * Transform the given real function, sampled on the given interval.
069         * <p>
070         * The formula is F<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; nk/N)
071         * </p>
072         *
073         * @param f the function to be sampled and transformed
074         * @param min the lower bound for the interval
075         * @param max the upper bound for the interval
076         * @param n the number of sample points
077         * @return the real transformed array
078         * @throws FunctionEvaluationException if function cannot be evaluated
079         * at some point
080         * @throws IllegalArgumentException if any parameters are invalid
081         */
082        public double[] transform(UnivariateRealFunction f,
083                                  double min, double max, int n)
084            throws FunctionEvaluationException, IllegalArgumentException {
085    
086            double data[] = FastFourierTransformer.sample(f, min, max, n);
087            data[0] = 0.0;
088            return fst(data);
089        }
090    
091        /**
092         * Transform the given real data set.
093         * <p>
094         * 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)
095         * </p>
096         *
097         * @param f the real data array to be transformed
098         * @return the real transformed array
099         * @throws IllegalArgumentException if any parameters are invalid
100         */
101        public double[] transform2(double f[]) throws IllegalArgumentException {
102    
103            double scaling_coefficient = FastMath.sqrt(2.0 / f.length);
104            return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
105        }
106    
107        /**
108         * Transform the given real function, sampled on the given interval.
109         * <p>
110         * 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)
111         * </p>
112         *
113         * @param f the function to be sampled and transformed
114         * @param min the lower bound for the interval
115         * @param max the upper bound for the interval
116         * @param n the number of sample points
117         * @return the real transformed array
118         * @throws FunctionEvaluationException if function cannot be evaluated
119         * at some point
120         * @throws IllegalArgumentException if any parameters are invalid
121         */
122        public double[] transform2(
123            UnivariateRealFunction f, double min, double max, int n)
124            throws FunctionEvaluationException, IllegalArgumentException {
125    
126            double data[] = FastFourierTransformer.sample(f, min, max, n);
127            data[0] = 0.0;
128            double scaling_coefficient = FastMath.sqrt(2.0 / n);
129            return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
130        }
131    
132        /**
133         * Inversely transform the given real data set.
134         * <p>
135         * 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)
136         * </p>
137         *
138         * @param f the real data array to be inversely transformed
139         * @return the real inversely transformed array
140         * @throws IllegalArgumentException if any parameters are invalid
141         */
142        public double[] inversetransform(double f[]) throws IllegalArgumentException {
143    
144            double scaling_coefficient = 2.0 / f.length;
145            return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient);
146        }
147    
148        /**
149         * Inversely transform the given real function, sampled on the given interval.
150         * <p>
151         * 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)
152         * </p>
153         *
154         * @param f the function to be sampled and inversely transformed
155         * @param min the lower bound for the interval
156         * @param max the upper bound for the interval
157         * @param n the number of sample points
158         * @return the real inversely transformed array
159         * @throws FunctionEvaluationException if function cannot be evaluated at some point
160         * @throws IllegalArgumentException if any parameters are invalid
161         */
162        public double[] inversetransform(UnivariateRealFunction f, double min, double max, int n)
163            throws FunctionEvaluationException, IllegalArgumentException {
164    
165            double data[] = FastFourierTransformer.sample(f, min, max, n);
166            data[0] = 0.0;
167            double scaling_coefficient = 2.0 / n;
168            return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient);
169        }
170    
171        /**
172         * Inversely transform the given real data set.
173         * <p>
174         * 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)
175         * </p>
176         *
177         * @param f the real data array to be inversely transformed
178         * @return the real inversely transformed array
179         * @throws IllegalArgumentException if any parameters are invalid
180         */
181        public double[] inversetransform2(double f[]) throws IllegalArgumentException {
182    
183            return transform2(f);
184        }
185    
186        /**
187         * Inversely transform the given real function, sampled on the given interval.
188         * <p>
189         * 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)
190         * </p>
191         *
192         * @param f the function to be sampled and inversely transformed
193         * @param min the lower bound for the interval
194         * @param max the upper bound for the interval
195         * @param n the number of sample points
196         * @return the real inversely transformed array
197         * @throws FunctionEvaluationException if function cannot be evaluated 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                        LocalizedFormats.FIRST_ELEMENT_NOT_ZERO,
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 = FastMath.sin(i * FastMath.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    }