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 */ 017package org.apache.commons.math4.transform; 018 019import java.util.Arrays; 020import java.util.function.DoubleUnaryOperator; 021 022import org.apache.commons.numbers.core.ArithmeticUtils; 023import org.apache.commons.numbers.complex.Complex; 024 025/** 026 * Implements the Fast Fourier Transform for transformation of one-dimensional 027 * real or complex data sets. For reference, see <em>Applied Numerical Linear 028 * Algebra</em>, ISBN 0898713897, chapter 6. 029 * <p> 030 * There are several variants of the discrete Fourier transform, with various 031 * normalization conventions, which are specified by the parameter 032 * {@link Norm}. 033 * <p> 034 * The current implementation of the discrete Fourier transform as a fast 035 * Fourier transform requires the length of the data set to be a power of 2. 036 * This greatly simplifies and speeds up the code. Users can pad the data with 037 * zeros to meet this requirement. There are other flavors of FFT, for 038 * reference, see S. Winograd, 039 * <i>On computing the discrete Fourier transform</i>, Mathematics of 040 * Computation, 32 (1978), 175 - 199. 041 */ 042public class FastFourierTransform implements ComplexTransform { 043 /** Number of array slots: 1 for "real" parts 1 for "imaginary" parts. */ 044 private static final int NUM_PARTS = 2; 045 /** 046 * {@code W_SUB_N_R[i]} is the real part of 047 * {@code exp(- 2 * i * pi / n)}: 048 * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}. 049 */ 050 private static final double[] W_SUB_N_R = { 051 0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1, 052 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1, 053 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1, 054 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1, 055 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1, 056 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1, 057 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1, 058 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0, 059 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 060 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 061 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 062 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 063 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 064 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 065 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, 066 0x1.0p0, 0x1.0p0, 0x1.0p0 }; 067 068 /** 069 * {@code W_SUB_N_I[i]} is the imaginary part of 070 * {@code exp(- 2 * i * pi / n)}: 071 * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}. 072 */ 073 private static final double[] W_SUB_N_I = { 074 0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1, 075 -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5, 076 -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9, 077 -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13, 078 -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17, 079 -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21, 080 -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25, 081 -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29, 082 -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33, 083 -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37, 084 -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41, 085 -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45, 086 -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49, 087 -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53, 088 -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57, 089 -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 }; 090 091 /** Type of DFT. */ 092 private final Norm normalization; 093 /** Inverse or forward. */ 094 private final boolean inverse; 095 096 /** 097 * @param normalization Normalization to be applied to the 098 * transformed data. 099 * @param inverse Whether to perform the inverse transform. 100 */ 101 public FastFourierTransform(final Norm normalization, 102 final boolean inverse) { 103 this.normalization = normalization; 104 this.inverse = inverse; 105 } 106 107 /** 108 * @param normalization Normalization to be applied to the 109 * transformed data. 110 */ 111 public FastFourierTransform(final Norm normalization) { 112 this(normalization, false); 113 } 114 115 /** 116 * Computes the standard transform of the data. 117 * Computation is done in place. 118 * Assumed layout of the input data: 119 * <ul> 120 * <li>{@code dataRI[0][i]}: Real part of the {@code i}-th data point,</li> 121 * <li>{@code dataRI[1][i]}: Imaginary part of the {@code i}-th data point.</li> 122 * </ul> 123 * 124 * @param dataRI Two-dimensional array of real and imaginary parts of the data. 125 * @throws IllegalArgumentException if the number of data points is not 126 * a power of two, if the number of rows of the specified array is not two, 127 * or the array is not rectangular. 128 */ 129 public void transformInPlace(final double[][] dataRI) { 130 if (dataRI.length != NUM_PARTS) { 131 throw new TransformException(TransformException.SIZE_MISMATCH, 132 dataRI.length, NUM_PARTS); 133 } 134 final double[] dataR = dataRI[0]; 135 final double[] dataI = dataRI[1]; 136 if (dataR.length != dataI.length) { 137 throw new TransformException(TransformException.SIZE_MISMATCH, 138 dataI.length, dataR.length); 139 } 140 final int n = dataR.length; 141 if (!ArithmeticUtils.isPowerOfTwo(n)) { 142 throw new TransformException(TransformException.NOT_POWER_OF_TWO, 143 Integer.valueOf(n)); 144 } 145 146 if (n == 1) { 147 return; 148 } else if (n == 2) { 149 final double srcR0 = dataR[0]; 150 final double srcI0 = dataI[0]; 151 final double srcR1 = dataR[1]; 152 final double srcI1 = dataI[1]; 153 154 // X_0 = x_0 + x_1 155 dataR[0] = srcR0 + srcR1; 156 dataI[0] = srcI0 + srcI1; 157 // X_1 = x_0 - x_1 158 dataR[1] = srcR0 - srcR1; 159 dataI[1] = srcI0 - srcI1; 160 161 normalizeTransformedData(dataRI); 162 return; 163 } 164 165 bitReversalShuffle2(dataR, dataI); 166 167 // Do 4-term DFT. 168 if (inverse) { 169 for (int i0 = 0; i0 < n; i0 += 4) { 170 final int i1 = i0 + 1; 171 final int i2 = i0 + 2; 172 final int i3 = i0 + 3; 173 174 final double srcR0 = dataR[i0]; 175 final double srcI0 = dataI[i0]; 176 final double srcR1 = dataR[i2]; 177 final double srcI1 = dataI[i2]; 178 final double srcR2 = dataR[i1]; 179 final double srcI2 = dataI[i1]; 180 final double srcR3 = dataR[i3]; 181 final double srcI3 = dataI[i3]; 182 183 // 4-term DFT 184 // X_0 = x_0 + x_1 + x_2 + x_3 185 dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3; 186 dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3; 187 // X_1 = x_0 - x_2 + j * (x_3 - x_1) 188 dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1); 189 dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3); 190 // X_2 = x_0 - x_1 + x_2 - x_3 191 dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3; 192 dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3; 193 // X_3 = x_0 - x_2 + j * (x_1 - x_3) 194 dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3); 195 dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1); 196 } 197 } else { 198 for (int i0 = 0; i0 < n; i0 += 4) { 199 final int i1 = i0 + 1; 200 final int i2 = i0 + 2; 201 final int i3 = i0 + 3; 202 203 final double srcR0 = dataR[i0]; 204 final double srcI0 = dataI[i0]; 205 final double srcR1 = dataR[i2]; 206 final double srcI1 = dataI[i2]; 207 final double srcR2 = dataR[i1]; 208 final double srcI2 = dataI[i1]; 209 final double srcR3 = dataR[i3]; 210 final double srcI3 = dataI[i3]; 211 212 // 4-term DFT 213 // X_0 = x_0 + x_1 + x_2 + x_3 214 dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3; 215 dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3; 216 // X_1 = x_0 - x_2 + j * (x_3 - x_1) 217 dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3); 218 dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1); 219 // X_2 = x_0 - x_1 + x_2 - x_3 220 dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3; 221 dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3; 222 // X_3 = x_0 - x_2 + j * (x_1 - x_3) 223 dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1); 224 dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3); 225 } 226 } 227 228 int lastN0 = 4; 229 int lastLogN0 = 2; 230 while (lastN0 < n) { 231 final int n0 = lastN0 << 1; 232 final int logN0 = lastLogN0 + 1; 233 final double wSubN0R = W_SUB_N_R[logN0]; 234 double wSubN0I = W_SUB_N_I[logN0]; 235 if (inverse) { 236 wSubN0I = -wSubN0I; 237 } 238 239 // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2). 240 for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) { 241 final int destOddStartIndex = destEvenStartIndex + lastN0; 242 243 double wSubN0ToRR = 1; 244 double wSubN0ToRI = 0; 245 246 for (int r = 0; r < lastN0; r++) { 247 final int destEvenStartIndexPlusR = destEvenStartIndex + r; 248 final int destOddStartIndexPlusR = destOddStartIndex + r; 249 250 final double grR = dataR[destEvenStartIndexPlusR]; 251 final double grI = dataI[destEvenStartIndexPlusR]; 252 final double hrR = dataR[destOddStartIndexPlusR]; 253 final double hrI = dataI[destOddStartIndexPlusR]; 254 255 final double a = wSubN0ToRR * hrR - wSubN0ToRI * hrI; 256 final double b = wSubN0ToRR * hrI + wSubN0ToRI * hrR; 257 // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr 258 dataR[destEvenStartIndexPlusR] = grR + a; 259 dataI[destEvenStartIndexPlusR] = grI + b; 260 // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr 261 dataR[destOddStartIndexPlusR] = grR - a; 262 dataI[destOddStartIndexPlusR] = grI - b; 263 264 // WsubN0ToR *= WsubN0R 265 final double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I; 266 final double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R; 267 wSubN0ToRR = nextWsubN0ToRR; 268 wSubN0ToRI = nextWsubN0ToRI; 269 } 270 } 271 272 lastN0 = n0; 273 lastLogN0 = logN0; 274 } 275 276 normalizeTransformedData(dataRI); 277 } 278 279 /** 280 * {@inheritDoc} 281 * 282 * @throws IllegalArgumentException if the length of the data array is not a power of two. 283 */ 284 @Override 285 public Complex[] apply(final double[] f) { 286 final double[][] dataRI = { 287 Arrays.copyOf(f, f.length), 288 new double[f.length] 289 }; 290 transformInPlace(dataRI); 291 return TransformUtils.createComplex(dataRI); 292 } 293 294 /** 295 * {@inheritDoc} 296 * 297 * @throws IllegalArgumentException if the number of sample points 298 * {@code n} is not a power of two, if the lower bound is greater than, 299 * or equal to the upper bound, if the number of sample points {@code n} 300 * is negative 301 */ 302 @Override 303 public Complex[] apply(final DoubleUnaryOperator f, 304 final double min, 305 final double max, 306 final int n) { 307 return apply(TransformUtils.sample(f, min, max, n)); 308 } 309 310 /** 311 * {@inheritDoc} 312 * 313 * @throws IllegalArgumentException if the length of the data array is 314 * not a power of two. 315 */ 316 @Override 317 public Complex[] apply(final Complex[] f) { 318 final double[][] dataRI = TransformUtils.createRealImaginary(f); 319 transformInPlace(dataRI); 320 return TransformUtils.createComplex(dataRI); 321 } 322 323 /** 324 * Applies normalization to the transformed data. 325 * 326 * @param dataRI Unscaled transformed data. 327 */ 328 private void normalizeTransformedData(final double[][] dataRI) { 329 final double[] dataR = dataRI[0]; 330 final double[] dataI = dataRI[1]; 331 final int n = dataR.length; 332 333 switch (normalization) { 334 case STD: 335 if (inverse) { 336 final double scaleFactor = 1d / n; 337 for (int i = 0; i < n; i++) { 338 dataR[i] *= scaleFactor; 339 dataI[i] *= scaleFactor; 340 } 341 } 342 343 break; 344 345 case UNIT: 346 final double scaleFactor = 1d / Math.sqrt(n); 347 for (int i = 0; i < n; i++) { 348 dataR[i] *= scaleFactor; 349 dataI[i] *= scaleFactor; 350 } 351 352 break; 353 354 default: 355 throw new IllegalStateException(); // Should never happen. 356 } 357 } 358 359 /** 360 * Performs identical index bit reversal shuffles on two arrays of 361 * identical size. 362 * Each element in the array is swapped with another element based 363 * on the bit-reversal of the index. 364 * For example, in an array with length 16, item at binary index 0011 365 * (decimal 3) would be swapped with the item at binary index 1100 366 * (decimal 12). 367 * 368 * @param a Array to be shuffled. 369 * @param b Array to be shuffled. 370 */ 371 private static void bitReversalShuffle2(double[] a, 372 double[] b) { 373 final int n = a.length; 374 final int halfOfN = n >> 1; 375 376 int j = 0; 377 for (int i = 0; i < n; i++) { 378 if (i < j) { 379 // swap indices i & j 380 double temp = a[i]; 381 a[i] = a[j]; 382 a[j] = temp; 383 384 temp = b[i]; 385 b[i] = b[j]; 386 b[j] = temp; 387 } 388 389 int k = halfOfN; 390 while (k <= j && k > 0) { 391 j -= k; 392 k >>= 1; 393 } 394 j += k; 395 } 396 } 397 398 /** 399 * Normalization types. 400 */ 401 public enum Norm { 402 /** 403 * Should be passed to the constructor of {@link FastFourierTransform} 404 * to use the <em>standard</em> normalization convention. This normalization 405 * convention is defined as follows 406 * <ul> 407 * <li>forward transform: \( y_n = \sum_{k = 0}^{N - 1} x_k e^{-2 \pi i n k / N} \),</li> 408 * <li>inverse transform: \( x_k = \frac{1}{N} \sum_{n = 0}^{N - 1} y_n e^{2 \pi i n k / N} \),</li> 409 * </ul> 410 * where \( N \) is the size of the data sample. 411 */ 412 STD, 413 414 /** 415 * Should be passed to the constructor of {@link FastFourierTransform} 416 * to use the <em>unitary</em> normalization convention. This normalization 417 * convention is defined as follows 418 * <ul> 419 * <li>forward transform: \( y_n = \frac{1}{\sqrt{N}} \sum_{k = 0}^{N - 1} x_k e^{-2 \pi i n k / N} \),</li> 420 * <li>inverse transform: \( x_k = \frac{1}{\sqrt{N}} \sum_{n = 0}^{N - 1} y_n e^{2 \pi i n k / N} \),</li> 421 * </ul> 422 * where \( N \) is the size of the data sample. 423 */ 424 UNIT; 425 } 426}