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 Cosine 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 * FCT 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 * Different from FFT and FST, FCT requires the length of data set to be 036 * power of 2 plus one. Users should especially pay attention to the 037 * function transformation on how this affects the sampling.</p> 038 * <p>As of version 2.0 this no longer implements Serializable</p> 039 * 040 * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $ 041 * @since 1.2 042 */ 043 public class FastCosineTransformer implements RealTransformer { 044 045 /** 046 * Construct a default transformer. 047 */ 048 public FastCosineTransformer() { 049 super(); 050 } 051 052 /** 053 * Transform the given real data set. 054 * <p> 055 * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 056 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 057 * </p> 058 * 059 * @param f the real data array to be transformed 060 * @return the real transformed array 061 * @throws IllegalArgumentException if any parameters are invalid 062 */ 063 public double[] transform(double f[]) throws IllegalArgumentException { 064 return fct(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> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 071 * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 072 * </p> 073 * 074 * @param f the function to be sampled and transformed 075 * @param min the lower bound for the interval 076 * @param max the upper bound for the interval 077 * @param n the number of sample points 078 * @return the real transformed array 079 * @throws FunctionEvaluationException if function cannot be evaluated 080 * at some point 081 * @throws IllegalArgumentException if any parameters are invalid 082 */ 083 public double[] transform(UnivariateRealFunction f, 084 double min, double max, int n) 085 throws FunctionEvaluationException, IllegalArgumentException { 086 double data[] = FastFourierTransformer.sample(f, min, max, n); 087 return fct(data); 088 } 089 090 /** 091 * Transform the given real data set. 092 * <p> 093 * The formula is F<sub>n</sub> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 094 * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π 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-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 }