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 */ 17package org.apache.commons.math.transform; 18 19import org.apache.commons.math.FunctionEvaluationException; 20import org.apache.commons.math.MathRuntimeException; 21import org.apache.commons.math.analysis.UnivariateRealFunction; 22import org.apache.commons.math.complex.Complex; 23import org.apache.commons.math.exception.util.LocalizedFormats; 24import org.apache.commons.math.util.FastMath; 25 26/** 27 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/ 28 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine Transform</a> 29 * for transformation of one-dimensional data sets. For reference, see 30 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3. 31 * <p> 32 * FCT is its own inverse, up to a multiplier depending on conventions. 33 * The equations are listed in the comments of the corresponding methods.</p> 34 * <p> 35 * Different from FFT and FST, FCT requires the length of data set to be 36 * power of 2 plus one. Users should especially pay attention to the 37 * function transformation on how this affects the sampling.</p> 38 * <p>As of version 2.0 this no longer implements Serializable</p> 39 * 40 * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $ 41 * @since 1.2 42 */ 43public class FastCosineTransformer implements RealTransformer { 44 45 /** 46 * Construct a default transformer. 47 */ 48 public FastCosineTransformer() { 49 super(); 50 } 51 52 /** 53 * Transform the given real data set. 54 * <p> 55 * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 56 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 57 * </p> 58 * 59 * @param f the real data array to be transformed 60 * @return the real transformed array 61 * @throws IllegalArgumentException if any parameters are invalid 62 */ 63 public double[] transform(double f[]) throws IllegalArgumentException { 64 return fct(f); 65 } 66 67 /** 68 * Transform the given real function, sampled on the given interval. 69 * <p> 70 * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 71 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 72 * </p> 73 * 74 * @param f the function to be sampled and transformed 75 * @param min the lower bound for the interval 76 * @param max the upper bound for the interval 77 * @param n the number of sample points 78 * @return the real transformed array 79 * @throws FunctionEvaluationException if function cannot be evaluated 80 * at some point 81 * @throws IllegalArgumentException if any parameters are invalid 82 */ 83 public double[] transform(UnivariateRealFunction f, 84 double min, double max, int n) 85 throws FunctionEvaluationException, IllegalArgumentException { 86 double data[] = FastFourierTransformer.sample(f, min, max, n); 87 return fct(data); 88 } 89 90 /** 91 * Transform the given real data set. 92 * <p> 93 * The formula is F<sub>n</sub> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 94 * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 95 * </p> 96 * 97 * @param f the real data array to be transformed 98 * @return the real transformed array 99 * @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-1)); 104 return FastFourierTransformer.scaleArray(fct(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> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 111 * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 112 * 113 * </p> 114 * 115 * @param f the function to be sampled and transformed 116 * @param min the lower bound for the interval 117 * @param max the upper bound for the interval 118 * @param n the number of sample points 119 * @return the real transformed array 120 * @throws FunctionEvaluationException if function cannot be evaluated 121 * at some point 122 * @throws IllegalArgumentException if any parameters are invalid 123 */ 124 public double[] transform2(UnivariateRealFunction f, 125 double min, double max, int n) 126 throws FunctionEvaluationException, IllegalArgumentException { 127 128 double data[] = FastFourierTransformer.sample(f, min, max, n); 129 double scaling_coefficient = FastMath.sqrt(2.0 / (n-1)); 130 return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient); 131 } 132 133 /** 134 * Inversely transform the given real data set. 135 * <p> 136 * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 137 * (2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 138 * </p> 139 * 140 * @param f the real data array to be inversely transformed 141 * @return the real inversely transformed array 142 * @throws IllegalArgumentException if any parameters are invalid 143 */ 144 public double[] inversetransform(double f[]) throws IllegalArgumentException { 145 146 double scaling_coefficient = 2.0 / (f.length - 1); 147 return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient); 148 } 149 150 /** 151 * Inversely transform the given real function, sampled on the given interval. 152 * <p> 153 * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 154 * (2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 155 * </p> 156 * 157 * @param f the function to be sampled and inversely transformed 158 * @param min the lower bound for the interval 159 * @param max the upper bound for the interval 160 * @param n the number of sample points 161 * @return the real inversely transformed array 162 * @throws FunctionEvaluationException if function cannot be evaluated at some point 163 * @throws IllegalArgumentException if any parameters are invalid 164 */ 165 public double[] inversetransform(UnivariateRealFunction f, 166 double min, double max, int n) 167 throws FunctionEvaluationException, IllegalArgumentException { 168 169 double data[] = FastFourierTransformer.sample(f, min, max, n); 170 double scaling_coefficient = 2.0 / (n - 1); 171 return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient); 172 } 173 174 /** 175 * Inversely transform the given real data set. 176 * <p> 177 * The formula is f<sub>k</sub> = √(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 178 * √(2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 179 * </p> 180 * 181 * @param f the real data array to be inversely transformed 182 * @return the real inversely transformed array 183 * @throws IllegalArgumentException if any parameters are invalid 184 */ 185 public double[] inversetransform2(double f[]) throws IllegalArgumentException { 186 return transform2(f); 187 } 188 189 /** 190 * Inversely transform the given real function, sampled on the given interval. 191 * <p> 192 * The formula is f<sub>k</sub> = √(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 193 * √(2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 194 * </p> 195 * 196 * @param f the function to be sampled and inversely transformed 197 * @param min the lower bound for the interval 198 * @param max the upper bound for the interval 199 * @param n the number of sample points 200 * @return the real inversely transformed array 201 * @throws FunctionEvaluationException if function cannot be evaluated 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 LocalizedFormats.NOT_POWER_OF_TWO_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 = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n-i]); 243 final double c = FastMath.cos(i * FastMath.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} 263