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 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 }