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.math3.transform;
018
019 import java.io.Serializable;
020
021 import org.apache.commons.math3.analysis.FunctionUtils;
022 import org.apache.commons.math3.analysis.UnivariateFunction;
023 import org.apache.commons.math3.complex.Complex;
024 import org.apache.commons.math3.exception.MathIllegalArgumentException;
025 import org.apache.commons.math3.exception.util.LocalizedFormats;
026 import org.apache.commons.math3.util.ArithmeticUtils;
027 import org.apache.commons.math3.util.FastMath;
028
029 /**
030 * Implements the Fast Cosine Transform for transformation of one-dimensional
031 * real data sets. For reference, see James S. Walker, <em>Fast Fourier
032 * Transforms</em>, chapter 3 (ISBN 0849371635).
033 * <p>
034 * There are several variants of the discrete cosine transform. The present
035 * implementation corresponds to DCT-I, with various normalization conventions,
036 * which are specified by the parameter {@link DctNormalization}.
037 * <p>
038 * DCT-I is equivalent to DFT of an <em>even extension</em> of the data series.
039 * More precisely, if x<sub>0</sub>, …, x<sub>N-1</sub> is the data set
040 * to be cosine transformed, the extended data set
041 * x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup>
042 * is defined as follows
043 * <ul>
044 * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 0 ≤ k < N,</li>
045 * <li>x<sub>k</sub><sup>#</sup> = x<sub>2N-2-k</sub>
046 * if N ≤ k < 2N - 2.</li>
047 * </ul>
048 * <p>
049 * Then, the standard DCT-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real
050 * data set x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em>
051 * of the N first elements of the DFT of the extended data set
052 * x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup>
053 * <br/>
054 * y<sub>n</sub> = (1 / 2) ∑<sub>k=0</sub><sup>2N-3</sup>
055 * x<sub>k</sub><sup>#</sup> exp[-2πi nk / (2N - 2)]
056 * k = 0, …, N-1.
057 * <p>
058 * The present implementation of the discrete cosine transform as a fast cosine
059 * transform requires the length of the data set to be a power of two plus one
060 * (N = 2<sup>n</sup> + 1). Besides, it implicitly assumes
061 * that the sampled function is even.
062 *
063 * @version $Id: FastCosineTransformer.java 1385310 2012-09-16 16:32:10Z tn $
064 * @since 1.2
065 */
066 public class FastCosineTransformer implements RealTransformer, Serializable {
067
068 /** Serializable version identifier. */
069 static final long serialVersionUID = 20120212L;
070
071 /** The type of DCT to be performed. */
072 private final DctNormalization normalization;
073
074 /**
075 * Creates a new instance of this class, with various normalization
076 * conventions.
077 *
078 * @param normalization the type of normalization to be applied to the
079 * transformed data
080 */
081 public FastCosineTransformer(final DctNormalization normalization) {
082 this.normalization = normalization;
083 }
084
085 /**
086 * {@inheritDoc}
087 *
088 * @throws MathIllegalArgumentException if the length of the data array is
089 * not a power of two plus one
090 */
091 public double[] transform(final double[] f, final TransformType type)
092 throws MathIllegalArgumentException {
093 if (type == TransformType.FORWARD) {
094 if (normalization == DctNormalization.ORTHOGONAL_DCT_I) {
095 final double s = FastMath.sqrt(2.0 / (f.length - 1));
096 return TransformUtils.scaleArray(fct(f), s);
097 }
098 return fct(f);
099 }
100 final double s2 = 2.0 / (f.length - 1);
101 final double s1;
102 if (normalization == DctNormalization.ORTHOGONAL_DCT_I) {
103 s1 = FastMath.sqrt(s2);
104 } else {
105 s1 = s2;
106 }
107 return TransformUtils.scaleArray(fct(f), s1);
108 }
109
110 /**
111 * {@inheritDoc}
112 *
113 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
114 * if the lower bound is greater than, or equal to the upper bound
115 * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
116 * if the number of sample points is negative
117 * @throws MathIllegalArgumentException if the number of sample points is
118 * not a power of two plus one
119 */
120 public double[] transform(final UnivariateFunction f,
121 final double min, final double max, final int n,
122 final TransformType type) throws MathIllegalArgumentException {
123
124 final double[] data = FunctionUtils.sample(f, min, max, n);
125 return transform(data, type);
126 }
127
128 /**
129 * Perform the FCT algorithm (including inverse).
130 *
131 * @param f the real data array to be transformed
132 * @return the real transformed array
133 * @throws MathIllegalArgumentException if the length of the data array is
134 * not a power of two plus one
135 */
136 protected double[] fct(double[] f)
137 throws MathIllegalArgumentException {
138
139 final double[] transformed = new double[f.length];
140
141 final int n = f.length - 1;
142 if (!ArithmeticUtils.isPowerOfTwo(n)) {
143 throw new MathIllegalArgumentException(
144 LocalizedFormats.NOT_POWER_OF_TWO_PLUS_ONE,
145 Integer.valueOf(f.length));
146 }
147 if (n == 1) { // trivial case
148 transformed[0] = 0.5 * (f[0] + f[1]);
149 transformed[1] = 0.5 * (f[0] - f[1]);
150 return transformed;
151 }
152
153 // construct a new array and perform FFT on it
154 final double[] x = new double[n];
155 x[0] = 0.5 * (f[0] + f[n]);
156 x[n >> 1] = f[n >> 1];
157 // temporary variable for transformed[1]
158 double t1 = 0.5 * (f[0] - f[n]);
159 for (int i = 1; i < (n >> 1); i++) {
160 final double a = 0.5 * (f[i] + f[n - i]);
161 final double b = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n - i]);
162 final double c = FastMath.cos(i * FastMath.PI / n) * (f[i] - f[n - i]);
163 x[i] = a - b;
164 x[n - i] = a + b;
165 t1 += c;
166 }
167 FastFourierTransformer transformer;
168 transformer = new FastFourierTransformer(DftNormalization.STANDARD);
169 Complex[] y = transformer.transform(x, TransformType.FORWARD);
170
171 // reconstruct the FCT result for the original array
172 transformed[0] = y[0].getReal();
173 transformed[1] = t1;
174 for (int i = 1; i < (n >> 1); i++) {
175 transformed[2 * i] = y[i].getReal();
176 transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary();
177 }
178 transformed[n] = y[n >> 1].getReal();
179
180 return transformed;
181 }
182 }