1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.transform; 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.FunctionEvaluationException; 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MathRuntimeException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.UnivariateRealFunction; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.complex.Complex; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/ 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Cosine Transform</a> 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * for transformation of one-dimensional data sets. For reference, see 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3. 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * FCT is its own inverse, up to a multiplier depending on conventions. 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The equations are listed in the comments of the corresponding methods.</p> 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Different from FFT and FST, FCT requires the length of data set to be 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * power of 2 plus one. Users should especially pay attention to the 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * function transformation on how this affects the sampling.</p> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>As of version 2.0 this no longer implements Serializable</p> 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class FastCosineTransformer implements RealTransformer { 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Construct a default transformer. 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public FastCosineTransformer() { 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(); 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Transform the given real data set. 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the real data array to be transformed 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real transformed array 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] transform(double f[]) throws IllegalArgumentException { 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return fct(f); 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Transform the given real function, sampled on the given interval. 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is F<sub>n</sub> = (1/2) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the function to be sampled and transformed 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param min the lower bound for the interval 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param max the upper bound for the interval 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param n the number of sample points 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real transformed array 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws FunctionEvaluationException if function cannot be evaluated 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * at some point 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] transform(UnivariateRealFunction f, 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double min, double max, int n) 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws FunctionEvaluationException, IllegalArgumentException { 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double data[] = FastFourierTransformer.sample(f, min, max, n); 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return fct(data); 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Transform the given real data set. 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is F<sub>n</sub> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the real data array to be transformed 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real transformed array 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] transform2(double f[]) throws IllegalArgumentException { 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double scaling_coefficient = FastMath.sqrt(2.0 / (f.length-1)); 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient); 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Transform the given real function, sampled on the given interval. 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is F<sub>n</sub> = √(1/2N) [f<sub>0</sub> + (-1)<sup>n</sup> f<sub>N</sub>] + 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * √(2/N) ∑<sub>k=1</sub><sup>N-1</sup> f<sub>k</sub> cos(π nk/N) 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the function to be sampled and transformed 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param min the lower bound for the interval 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param max the upper bound for the interval 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param n the number of sample points 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real transformed array 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws FunctionEvaluationException if function cannot be evaluated 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * at some point 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] transform2(UnivariateRealFunction f, 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double min, double max, int n) 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws FunctionEvaluationException, IllegalArgumentException { 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double data[] = FastFourierTransformer.sample(f, min, max, n); 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double scaling_coefficient = FastMath.sqrt(2.0 / (n-1)); 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient); 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Inversely transform the given real data set. 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the real data array to be inversely transformed 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real inversely transformed array 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] inversetransform(double f[]) throws IllegalArgumentException { 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double scaling_coefficient = 2.0 / (f.length - 1); 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return FastFourierTransformer.scaleArray(fct(f), scaling_coefficient); 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Inversely transform the given real function, sampled on the given interval. 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is f<sub>k</sub> = (1/N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the function to be sampled and inversely transformed 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param min the lower bound for the interval 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param max the upper bound for the interval 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param n the number of sample points 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real inversely transformed array 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws FunctionEvaluationException if function cannot be evaluated at some point 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] inversetransform(UnivariateRealFunction f, 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double min, double max, int n) 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws FunctionEvaluationException, IllegalArgumentException { 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double data[] = FastFourierTransformer.sample(f, min, max, n); 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double scaling_coefficient = 2.0 / (n - 1); 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return FastFourierTransformer.scaleArray(fct(data), scaling_coefficient); 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Inversely transform the given real data set. 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is f<sub>k</sub> = √(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * √(2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the real data array to be inversely transformed 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real inversely transformed array 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] inversetransform2(double f[]) throws IllegalArgumentException { 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return transform2(f); 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Inversely transform the given real function, sampled on the given interval. 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The formula is f<sub>k</sub> = √(1/2N) [F<sub>0</sub> + (-1)<sup>k</sup> F<sub>N</sub>] + 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * √(2/N) ∑<sub>n=1</sub><sup>N-1</sup> F<sub>n</sub> cos(π nk/N) 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the function to be sampled and inversely transformed 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param min the lower bound for the interval 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param max the upper bound for the interval 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param n the number of sample points 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real inversely transformed array 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws FunctionEvaluationException if function cannot be evaluated at some point 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] inversetransform2(UnivariateRealFunction f, 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double min, double max, int n) 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws FunctionEvaluationException, IllegalArgumentException { 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return transform2(f, min, max, n); 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Perform the FCT algorithm (including inverse). 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f the real data array to be transformed 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the real transformed array 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if any parameters are invalid 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double[] fct(double f[]) 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException { 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double transformed[] = new double[f.length]; 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = f.length - 1; 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!FastFourierTransformer.isPowerOf2(n)) { 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.NOT_POWER_OF_TWO_PLUS_ONE, 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond f.length); 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (n == 1) { // trivial case 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[0] = 0.5 * (f[0] + f[1]); 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[1] = 0.5 * (f[0] - f[1]); 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return transformed; 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // construct a new array and perform FFT on it 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] x = new double[n]; 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond x[0] = 0.5 * (f[0] + f[n]); 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond x[n >> 1] = f[n >> 1]; 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double t1 = 0.5 * (f[0] - f[n]); // temporary variable for transformed[1] 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 1; i < (n >> 1); i++) { 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double a = 0.5 * (f[i] + f[n-i]); 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double b = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n-i]); 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double c = FastMath.cos(i * FastMath.PI / n) * (f[i] - f[n-i]); 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond x[i] = a - b; 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond x[n-i] = a + b; 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond t1 += c; 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastFourierTransformer transformer = new FastFourierTransformer(); 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Complex y[] = transformer.transform(x); 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // reconstruct the FCT result for the original array 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[0] = y[0].getReal(); 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[1] = t1; 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 1; i < (n >> 1); i++) { 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[2 * i] = y[i].getReal(); 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary(); 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond transformed[n] = y[n >> 1].getReal(); 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return transformed; 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 263