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.FunctionEvaluationException;
20 import org.apache.commons.math.MathRuntimeException;
21 import org.apache.commons.math.analysis.UnivariateRealFunction;
22 import org.apache.commons.math.complex.Complex;
23
24 /**
25 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/
26 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine 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 * FCT 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 * Different from FFT and FST, FCT requires the length of data set to be
34 * power of 2 plus one. Users should especially pay attention to the
35 * function transformation on how this affects the sampling.</p>
36 * <p>As of version 2.0 this no longer implements Serializable</p>
37 *
38 * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
39 * @since 1.2
40 */
41 public class FastCosineTransformer implements RealTransformer {
42
43 /**
44 * Construct a default transformer.
45 */
46 public FastCosineTransformer() {
47 super();
48 }
49
50 /**
51 * Transform the given real data set.
52 * <p>
53 * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
54 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N)
55 * </p>
56 *
57 * @param f the real data array to be transformed
58 * @return the real transformed array
59 * @throws IllegalArgumentException if any parameters are invalid
60 */
61 public double[] transform(double f[]) throws IllegalArgumentException {
62 return fct(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> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
69 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N)
70 * </p>
71 *
72 * @param f the function to be sampled and transformed
73 * @param min the lower bound for the interval
74 * @param max the upper bound for the interval
75 * @param n the number of sample points
76 * @return the real transformed array
77 * @throws FunctionEvaluationException if function cannot be evaluated
78 * at some point
79 * @throws IllegalArgumentException if any parameters are invalid
80 */
81 public double[] transform(UnivariateRealFunction f,
82 double min, double max, int n)
83 throws FunctionEvaluationException, IllegalArgumentException {
84 double data[] = FastFourierTransformer.sample(f, min, max, n);
85 return fct(data);
86 }
87
88 /**
89 * Transform the given real data set.
90 * <p>
91 * The formula is F<sub>n</sub> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] +
92 * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π 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-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 double A, B, C, F1, x[], F[] = new double[f.length];
222
223 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 F[0] = 0.5 * (f[0] + f[1]);
231 F[1] = 0.5 * (f[0] - f[1]);
232 return F;
233 }
234
235 // construct a new array and perform FFT on it
236 x = new double[N];
237 x[0] = 0.5 * (f[0] + f[N]);
238 x[N >> 1] = f[N >> 1];
239 F1 = 0.5 * (f[0] - f[N]); // temporary variable for F[1]
240 for (int i = 1; i < (N >> 1); i++) {
241 A = 0.5 * (f[i] + f[N-i]);
242 B = Math.sin(i * Math.PI / N) * (f[i] - f[N-i]);
243 C = Math.cos(i * Math.PI / N) * (f[i] - f[N-i]);
244 x[i] = A - B;
245 x[N-i] = A + B;
246 F1 += C;
247 }
248 FastFourierTransformer transformer = new FastFourierTransformer();
249 Complex y[] = transformer.transform(x);
250
251 // reconstruct the FCT result for the original array
252 F[0] = y[0].getReal();
253 F[1] = F1;
254 for (int i = 1; i < (N >> 1); i++) {
255 F[2*i] = y[i].getReal();
256 F[2*i+1] = F[2*i-1] - y[i].getImaginary();
257 }
258 F[N] = y[N >> 1].getReal();
259
260 return F;
261 }
262 }