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 Cosine 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 * FCT 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 * Different from FFT and FST, FCT requires the length of data set to be
034 * power of 2 plus one. Users should especially pay attention to the
035 * function transformation on how this affects the sampling.</p>
036 * <p>As of version 2.0 this no longer implements Serializable</p>
037 *
038 * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
039 * @since 1.2
040 */
041 public class FastCosineTransformer implements RealTransformer {
042
043 /**
044 * Construct a default transformer.
045 */
046 public FastCosineTransformer() {
047 super();
048 }
049
050 /**
051 * Transform the given real data set.
052 * <p>
053 * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
054 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N)
055 * </p>
056 *
057 * @param f the real data array to be transformed
058 * @return the real transformed array
059 * @throws IllegalArgumentException if any parameters are invalid
060 */
061 public double[] transform(double f[]) throws IllegalArgumentException {
062 return fct(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> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
069 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N)
070 * </p>
071 *
072 * @param f the function to be sampled and transformed
073 * @param min the lower bound for the interval
074 * @param max the upper bound for the interval
075 * @param n the number of sample points
076 * @return the real transformed array
077 * @throws FunctionEvaluationException if function cannot be evaluated
078 * at some point
079 * @throws IllegalArgumentException if any parameters are invalid
080 */
081 public double[] transform(UnivariateRealFunction f,
082 double min, double max, int n)
083 throws FunctionEvaluationException, IllegalArgumentException {
084 double data[] = FastFourierTransformer.sample(f, min, max, n);
085 return fct(data);
086 }
087
088 /**
089 * Transform the given real data set.
090 * <p>
091 * The formula is F<sub>n</sub> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
092 * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π 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-1));
102 return FastFourierTransformer.scaleArray(fct(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> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
109 * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N)
110 *
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(UnivariateRealFunction f,
123 double min, double max, int n)
124 throws FunctionEvaluationException, IllegalArgumentException {
125
126 double data[] = FastFourierTransformer.sample(f, min, max, n);
127 double scaling_coefficient = Math.sqrt(2.0 / (n-1));
128 return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
129 }
130
131 /**
132 * Inversely transform the given real data set.
133 * <p>
134 * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
135 * (2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π 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 - 1);
145 return FastFourierTransformer.scaleArray(fct(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> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
152 * (2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N)
153 * </p>
154 *
155 * @param f the function to be sampled and inversely transformed
156 * @param min the lower bound for the interval
157 * @param max the upper bound for the interval
158 * @param n the number of sample points
159 * @return the real inversely transformed array
160 * @throws FunctionEvaluationException if function cannot be evaluated
161 * at some point
162 * @throws IllegalArgumentException if any parameters are invalid
163 */
164 public double[] inversetransform(UnivariateRealFunction f,
165 double min, double max, int n)
166 throws FunctionEvaluationException, IllegalArgumentException {
167
168 double data[] = FastFourierTransformer.sample(f, min, max, n);
169 double scaling_coefficient = 2.0 / (n - 1);
170 return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient);
171 }
172
173 /**
174 * Inversely transform the given real data set.
175 * <p>
176 * The formula is f<sub>k</sub> = √(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
177 * √(2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N)
178 * </p>
179 *
180 * @param f the real data array to be inversely transformed
181 * @return the real inversely transformed array
182 * @throws IllegalArgumentException if any parameters are invalid
183 */
184 public double[] inversetransform2(double f[]) throws IllegalArgumentException {
185 return transform2(f);
186 }
187
188 /**
189 * Inversely transform the given real function, sampled on the given interval.
190 * <p>
191 * The formula is f<sub>k</sub> = √(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] +
192 * √(2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N)
193 * </p>
194 *
195 * @param f the function to be sampled and inversely transformed
196 * @param min the lower bound for the interval
197 * @param max the upper bound for the interval
198 * @param n the number of sample points
199 * @return the real inversely transformed array
200 * @throws FunctionEvaluationException if function cannot be evaluated
201 * at some point
202 * @throws IllegalArgumentException if any parameters are invalid
203 */
204 public double[] inversetransform2(UnivariateRealFunction f,
205 double min, double max, int n)
206 throws FunctionEvaluationException, IllegalArgumentException {
207
208 return transform2(f, min, max, n);
209 }
210
211 /**
212 * Perform the FCT algorithm (including inverse).
213 *
214 * @param f the real data array to be transformed
215 * @return the real transformed array
216 * @throws IllegalArgumentException if any parameters are invalid
217 */
218 protected double[] fct(double f[])
219 throws IllegalArgumentException {
220
221 final double transformed[] = new double[f.length];
222
223 final int n = f.length - 1;
224 if (!FastFourierTransformer.isPowerOf2(n)) {
225 throw MathRuntimeException.createIllegalArgumentException(
226 "{0} is not a power of 2 plus one",
227 f.length);
228 }
229 if (n == 1) { // trivial case
230 transformed[0] = 0.5 * (f[0] + f[1]);
231 transformed[1] = 0.5 * (f[0] - f[1]);
232 return transformed;
233 }
234
235 // construct a new array and perform FFT on it
236 final double[] x = new double[n];
237 x[0] = 0.5 * (f[0] + f[n]);
238 x[n >> 1] = f[n >> 1];
239 double t1 = 0.5 * (f[0] - f[n]); // temporary variable for transformed[1]
240 for (int i = 1; i < (n >> 1); i++) {
241 final double a = 0.5 * (f[i] + f[n-i]);
242 final double b = Math.sin(i * Math.PI / n) * (f[i] - f[n-i]);
243 final double c = Math.cos(i * Math.PI / n) * (f[i] - f[n-i]);
244 x[i] = a - b;
245 x[n-i] = a + b;
246 t1 += c;
247 }
248 FastFourierTransformer transformer = new FastFourierTransformer();
249 Complex y[] = transformer.transform(x);
250
251 // reconstruct the FCT result for the original array
252 transformed[0] = y[0].getReal();
253 transformed[1] = t1;
254 for (int i = 1; i < (n >> 1); i++) {
255 transformed[2 * i] = y[i].getReal();
256 transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary();
257 }
258 transformed[n] = y[n >> 1].getReal();
259
260 return transformed;
261 }
262 }