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.analysis.*;
020    import org.apache.commons.math.complex.*;
021    import org.apache.commons.math.FunctionEvaluationException;
022    import org.apache.commons.math.MathRuntimeException;
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: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 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> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; 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> = &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; 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> = &radic;(2/N) &sum;<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(&pi; 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> = &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    }