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 */ 17 package org.apache.commons.math4.transform; 18 19 import java.util.function.UnaryOperator; 20 import java.util.function.DoubleUnaryOperator; 21 22 import org.apache.commons.numbers.complex.Complex; 23 import org.apache.commons.numbers.core.ArithmeticUtils; 24 25 /** 26 * Implements the Fast Cosine Transform for transformation of one-dimensional 27 * real data sets. For reference, see James S. Walker, <em>Fast Fourier 28 * Transforms</em>, chapter 3 (ISBN 0849371635). 29 * <p> 30 * There are several variants of the discrete cosine transform. The present 31 * implementation corresponds to DCT-I, with various normalization conventions, 32 * which are specified by the parameter {@link Norm}. 33 * <p> 34 * DCT-I is equivalent to DFT of an <em>even extension</em> of the data series. 35 * More precisely, if x<sub>0</sub>, …, x<sub>N-1</sub> is the data set 36 * to be cosine transformed, the extended data set 37 * x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup> 38 * is defined as follows 39 * <ul> 40 * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 0 ≤ k < N,</li> 41 * <li>x<sub>k</sub><sup>#</sup> = x<sub>2N-2-k</sub> 42 * if N ≤ k < 2N - 2.</li> 43 * </ul> 44 * <p> 45 * Then, the standard DCT-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real 46 * data set x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em> 47 * of the N first elements of the DFT of the extended data set 48 * x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup> 49 * <br> 50 * y<sub>n</sub> = (1 / 2) ∑<sub>k=0</sub><sup>2N-3</sup> 51 * x<sub>k</sub><sup>#</sup> exp[-2πi nk / (2N - 2)] 52 * k = 0, …, N-1. 53 * <p> 54 * The present implementation of the discrete cosine transform as a fast cosine 55 * transform requires the length of the data set to be a power of two plus one 56 * (N = 2<sup>n</sup> + 1). Besides, it implicitly assumes 57 * that the sampled function is even. 58 */ 59 public class FastCosineTransform implements RealTransform { 60 /** Operation to be performed. */ 61 private final UnaryOperator<double[]> op; 62 63 /** 64 * @param normalization Normalization to be applied to the 65 * transformed data. 66 * @param inverse Whether to perform the inverse transform. 67 */ 68 public FastCosineTransform(final Norm normalization, 69 final boolean inverse) { 70 op = create(normalization, inverse); 71 } 72 73 /** 74 * @param normalization Normalization to be applied to the 75 * transformed data. 76 */ 77 public FastCosineTransform(final Norm normalization) { 78 this(normalization, false); 79 } 80 81 /** 82 * {@inheritDoc} 83 * 84 * @throws IllegalArgumentException if the length of the data array is 85 * not a power of two plus one. 86 */ 87 @Override 88 public double[] apply(final double[] f) { 89 return op.apply(f); 90 } 91 92 /** 93 * {@inheritDoc} 94 * 95 * @throws IllegalArgumentException if the number of sample points is 96 * not a power of two plus one, if the lower bound is greater than or 97 * equal to the upper bound, if the number of sample points is negative. 98 */ 99 @Override 100 public double[] apply(final DoubleUnaryOperator f, 101 final double min, 102 final double max, 103 final int n) { 104 return apply(TransformUtils.sample(f, min, max, n)); 105 } 106 107 /** 108 * Perform the FCT algorithm (including inverse). 109 * 110 * @param f Data to be transformed. 111 * @return the transformed array. 112 * @throws IllegalArgumentException if the length of the data array is 113 * not a power of two plus one. 114 */ 115 private double[] fct(double[] f) { 116 final int n = f.length - 1; 117 if (!ArithmeticUtils.isPowerOfTwo(n)) { 118 throw new TransformException(TransformException.NOT_POWER_OF_TWO_PLUS_ONE, 119 Integer.valueOf(f.length)); 120 } 121 122 final double[] transformed = new double[f.length]; 123 124 if (n == 1) { // trivial case 125 transformed[0] = 0.5 * (f[0] + f[1]); 126 transformed[1] = 0.5 * (f[0] - f[1]); 127 return transformed; 128 } 129 130 // construct a new array and perform FFT on it 131 final double[] x = new double[n]; 132 x[0] = 0.5 * (f[0] + f[n]); 133 final int nShifted = n >> 1; 134 x[nShifted] = f[nShifted]; 135 // temporary variable for transformed[1] 136 double t1 = 0.5 * (f[0] - f[n]); 137 final double piOverN = Math.PI / n; 138 for (int i = 1; i < nShifted; i++) { 139 final int nMi = n - i; 140 final double fi = f[i]; 141 final double fnMi = f[nMi]; 142 final double a = 0.5 * (fi + fnMi); 143 final double arg = i * piOverN; 144 final double b = Math.sin(arg) * (fi - fnMi); 145 final double c = Math.cos(arg) * (fi - fnMi); 146 x[i] = a - b; 147 x[nMi] = a + b; 148 t1 += c; 149 } 150 final FastFourierTransform transformer = new FastFourierTransform(FastFourierTransform.Norm.STD, 151 false); 152 final Complex[] y = transformer.apply(x); 153 154 // reconstruct the FCT result for the original array 155 transformed[0] = y[0].getReal(); 156 transformed[1] = t1; 157 for (int i = 1; i < nShifted; i++) { 158 final int i2 = 2 * i; 159 transformed[i2] = y[i].getReal(); 160 transformed[i2 + 1] = transformed[i2 - 1] - y[i].getImaginary(); 161 } 162 transformed[n] = y[nShifted].getReal(); 163 164 return transformed; 165 } 166 167 /** 168 * Factory method. 169 * 170 * @param normalization Normalization to be applied to the 171 * transformed data. 172 * @param inverse Whether to perform the inverse transform. 173 * @return the transform operator. 174 */ 175 private UnaryOperator<double[]> create(final Norm normalization, 176 final boolean inverse) { 177 if (inverse) { 178 return normalization == Norm.ORTHO ? 179 f -> TransformUtils.scaleInPlace(fct(f), Math.sqrt(2d / (f.length - 1))) : 180 f -> TransformUtils.scaleInPlace(fct(f), 2d / (f.length - 1)); 181 } else { 182 return normalization == Norm.ORTHO ? 183 f -> TransformUtils.scaleInPlace(fct(f), Math.sqrt(2d / (f.length - 1))) : 184 f -> fct(f); 185 } 186 } 187 188 /** 189 * Normalization types. 190 */ 191 public enum Norm { 192 /** 193 * Should be passed to the constructor of {@link FastCosineTransform} 194 * to use the <em>standard</em> normalization convention. The standard 195 * DCT-I normalization convention is defined as follows 196 * <ul> 197 * <li>forward transform: 198 * y<sub>n</sub> = (1/2) [x<sub>0</sub> + (-1)<sup>n</sup>x<sub>N-1</sub>] 199 * + ∑<sub>k=1</sub><sup>N-2</sup> 200 * x<sub>k</sub> cos[π nk / (N - 1)],</li> 201 * <li>inverse transform: 202 * x<sub>k</sub> = [1 / (N - 1)] [y<sub>0</sub> 203 * + (-1)<sup>k</sup>y<sub>N-1</sub>] 204 * + [2 / (N - 1)] ∑<sub>n=1</sub><sup>N-2</sup> 205 * y<sub>n</sub> cos[π nk / (N - 1)],</li> 206 * </ul> 207 * where N is the size of the data sample. 208 */ 209 STD, 210 211 /** 212 * Should be passed to the constructor of {@link FastCosineTransform} 213 * to use the <em>orthogonal</em> normalization convention. The orthogonal 214 * DCT-I normalization convention is defined as follows 215 * <ul> 216 * <li>forward transform: 217 * y<sub>n</sub> = [2(N - 1)]<sup>-1/2</sup> [x<sub>0</sub> 218 * + (-1)<sup>n</sup>x<sub>N-1</sub>] 219 * + [2 / (N - 1)]<sup>1/2</sup> ∑<sub>k=1</sub><sup>N-2</sup> 220 * x<sub>k</sub> cos[π nk / (N - 1)],</li> 221 * <li>inverse transform: 222 * x<sub>k</sub> = [2(N - 1)]<sup>-1/2</sup> [y<sub>0</sub> 223 * + (-1)<sup>k</sup>y<sub>N-1</sub>] 224 * + [2 / (N - 1)]<sup>1/2</sup> ∑<sub>n=1</sub><sup>N-2</sup> 225 * y<sub>n</sub> cos[π nk / (N - 1)],</li> 226 * </ul> 227 * which makes the transform orthogonal. N is the size of the data sample. 228 */ 229 ORTHO; 230 } 231 }