/* * Copyright (c) 2007, Oracle and/or its affiliates. All rights reserved. * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. * * This code is free software; you can redistribute it and/or modify it * under the terms of the GNU General Public License version 2 only, as * published by the Free Software Foundation. Oracle designates this * particular file as subject to the "Classpath" exception as provided * by Oracle in the LICENSE file that accompanied this code. * * This code is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License * version 2 for more details (a copy is included in the LICENSE file that * accompanied this code). * * You should have received a copy of the GNU General Public License version * 2 along with this work; if not, write to the Free Software Foundation, * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. * * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA * or visit www.oracle.com if you need additional information or have any * questions. */package com.sun.media.sound; /** * Fast Fourier Transformer. * * @author Karl Helgason */public final class FFT { private double[] w; private int fftFrameSize; private int sign; private int[] bitm_array; private int fftFrameSize2; // Sign = -1 is FFT, 1 is IFFT (inverse FFT) // Data = Interlaced double array to be transformed. // The order is: real (sin), complex (cos) // Framesize must be power of 2 public FFT(int fftFrameSize, int sign) { w = computeTwiddleFactors(fftFrameSize, sign); this.fftFrameSize = fftFrameSize; this.sign = sign; fftFrameSize2 = fftFrameSize << 1; // Pre-process Bit-Reversal bitm_array = new int[fftFrameSize2]; for (int i = 2; i < fftFrameSize2; i += 2) { int j; int bitm; for (bitm = 2, j = 0; bitm < fftFrameSize2; bitm <<= 1) { if ((i & bitm) != 0) j++; j <<= 1; } bitm_array[i] = j; } } public void transform(double[] data) { bitreversal(data); calc(fftFrameSize, data, sign, w); } private final static double[] computeTwiddleFactors(int fftFrameSize, int sign) { int imax = (int) (Math.log(fftFrameSize) / Math.log(2.)); double[] warray = new double[(fftFrameSize - 1) * 4]; int w_index = 0; for (int i = 0, nstep = 2; i < imax; i++) { int jmax = nstep; nstep <<= 1; double wr = 1.0; double wi = 0.0; double arg = Math.PI / (jmax >> 1); double wfr = Math.cos(arg); double wfi = sign * Math.sin(arg); for (int j = 0; j < jmax; j += 2) { warray[w_index++] = wr; warray[w_index++] = wi; double tempr = wr; wr = tempr * wfr - wi * wfi; wi = tempr * wfi + wi * wfr; } } // PRECOMPUTATION of wwr1, wwi1 for factor 4 Decomposition (3 * complex // operators and 8 +/- complex operators) { w_index = 0; int w_index2 = warray.length >> 1; for (int i = 0, nstep = 2; i < (imax - 1); i++) { int jmax = nstep; nstep *= 2; int ii = w_index + jmax; for (int j = 0; j < jmax; j += 2) { double wr = warray[w_index++]; double wi = warray[w_index++]; double wr1 = warray[ii++]; double wi1 = warray[ii++]; warray[w_index2++] = wr * wr1 - wi * wi1; warray[w_index2++] = wr * wi1 + wi * wr1; } } } return warray; } private final static void calc(int fftFrameSize, double[] data, int sign, double[] w) { final int fftFrameSize2 = fftFrameSize << 1; int nstep = 2; if (nstep >= fftFrameSize2) return; int i = nstep - 2; if (sign == -1) calcF4F(fftFrameSize, data, i, nstep, w); else calcF4I(fftFrameSize, data, i, nstep, w); } private final static void calcF2E(int fftFrameSize, double[] data, int i, int nstep, double[] w) { int jmax = nstep; for (int n = 0; n < jmax; n += 2) { double wr = w[i++]; double wi = w[i++]; int m = n + jmax; double datam_r = data[m]; double datam_i = data[m + 1]; double datan_r = data[n]; double datan_i = data[n + 1]; double tempr = datam_r * wr - datam_i * wi; double tempi = datam_r * wi + datam_i * wr; data[m] = datan_r - tempr; data[m + 1] = datan_i - tempi; data[n] = datan_r + tempr; data[n + 1] = datan_i + tempi; } return; } // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/- // complex operators private final static void calcF4F(int fftFrameSize, double[] data, int i, int nstep, double[] w) { final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize; // Factor-4 Decomposition int w_len = w.length >> 1; while (nstep < fftFrameSize2) { if (nstep << 2 == fftFrameSize2) { // Goto Factor-4 Final Decomposition // calcF4E(data, i, nstep, -1, w); calcF4FE(fftFrameSize, data, i, nstep, w); return; } int jmax = nstep; int nnstep = nstep << 1; if (nnstep == fftFrameSize2) { // Factor-4 Decomposition not possible calcF2E(fftFrameSize, data, i, nstep, w); return; } nstep <<= 2; int ii = i + jmax; int iii = i + w_len; { i += 2; ii += 2; iii += 2; for (int n = 0; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r; double tempi = datam1_i; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r; double n2w1i = datan2_i; double m2ww1r = datam2_r; double m2ww1i = datam2_i; tempr = m2ww1r - n2w1r; tempi = m2ww1i - n2w1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } for (int j = 2; j < jmax; j += 2) { double wr = w[i++]; double wi = w[i++]; double wr1 = w[ii++]; double wi1 = w[ii++]; double wwr1 = w[iii++]; double wwi1 = w[iii++]; // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be // precomputed!!! // double wwi1 = wr * wi1 + wi * wr1; for (int n = j; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r * wr - datam1_i * wi; double tempi = datam1_r * wi + datam1_i * wr; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r * wr1 - datan2_i * wi1; double n2w1i = datan2_r * wi1 + datan2_i * wr1; double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1; double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1; tempr = m2ww1r - n2w1r; tempi = m2ww1i - n2w1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } i += jmax << 1; } calcF2E(fftFrameSize, data, i, nstep, w); } // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/- // complex operators private final static void calcF4I(int fftFrameSize, double[] data, int i, int nstep, double[] w) { final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize; // Factor-4 Decomposition int w_len = w.length >> 1; while (nstep < fftFrameSize2) { if (nstep << 2 == fftFrameSize2) { // Goto Factor-4 Final Decomposition // calcF4E(data, i, nstep, 1, w); calcF4IE(fftFrameSize, data, i, nstep, w); return; } int jmax = nstep; int nnstep = nstep << 1; if (nnstep == fftFrameSize2) { // Factor-4 Decomposition not possible calcF2E(fftFrameSize, data, i, nstep, w); return; } nstep <<= 2; int ii = i + jmax; int iii = i + w_len; { i += 2; ii += 2; iii += 2; for (int n = 0; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r; double tempi = datam1_i; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r; double n2w1i = datan2_i; double m2ww1r = datam2_r; double m2ww1i = datam2_i; tempr = n2w1r - m2ww1r; tempi = n2w1i - m2ww1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } for (int j = 2; j < jmax; j += 2) { double wr = w[i++]; double wi = w[i++]; double wr1 = w[ii++]; double wi1 = w[ii++]; double wwr1 = w[iii++]; double wwi1 = w[iii++]; // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be // precomputed!!! // double wwi1 = wr * wi1 + wi * wr1; for (int n = j; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r * wr - datam1_i * wi; double tempi = datam1_r * wi + datam1_i * wr; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r * wr1 - datan2_i * wi1; double n2w1i = datan2_r * wi1 + datan2_i * wr1; double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1; double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1; tempr = n2w1r - m2ww1r; tempi = n2w1i - m2ww1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } i += jmax << 1; } calcF2E(fftFrameSize, data, i, nstep, w); } // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-Lines 1–500 of 748. Load more
/* * Copyright (c) 2007, Oracle and/or its affiliates. All rights reserved. * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. * * This code is free software; you can redistribute it and/or modify it * under the terms of the GNU General Public License version 2 only, as * published by the Free Software Foundation. Oracle designates this * particular file as subject to the "Classpath" exception as provided * by Oracle in the LICENSE file that accompanied this code. * * This code is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License * version 2 for more details (a copy is included in the LICENSE file that * accompanied this code). * * You should have received a copy of the GNU General Public License version * 2 along with this work; if not, write to the Free Software Foundation, * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. * * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA * or visit www.oracle.com if you need additional information or have any * questions. */ package com.sun.media.sound; /** * Fast Fourier Transformer. * * @author Karl Helgason */ public final class FFT { private double[] w; private int fftFrameSize; private int sign; private int[] bitm_array; private int fftFrameSize2; // Sign = -1 is FFT, 1 is IFFT (inverse FFT) // Data = Interlaced double array to be transformed. // The order is: real (sin), complex (cos) // Framesize must be power of 2 public FFT(int fftFrameSize, int sign) { w = computeTwiddleFactors(fftFrameSize, sign); this.fftFrameSize = fftFrameSize; this.sign = sign; fftFrameSize2 = fftFrameSize << 1; // Pre-process Bit-Reversal bitm_array = new int[fftFrameSize2]; for (int i = 2; i < fftFrameSize2; i += 2) { int j; int bitm; for (bitm = 2, j = 0; bitm < fftFrameSize2; bitm <<= 1) { if ((i & bitm) != 0) j++; j <<= 1; } bitm_array[i] = j; } } public void transform(double[] data) { bitreversal(data); calc(fftFrameSize, data, sign, w); } private final static double[] computeTwiddleFactors(int fftFrameSize, int sign) { int imax = (int) (Math.log(fftFrameSize) / Math.log(2.)); double[] warray = new double[(fftFrameSize - 1) * 4]; int w_index = 0; for (int i = 0, nstep = 2; i < imax; i++) { int jmax = nstep; nstep <<= 1; double wr = 1.0; double wi = 0.0; double arg = Math.PI / (jmax >> 1); double wfr = Math.cos(arg); double wfi = sign * Math.sin(arg); for (int j = 0; j < jmax; j += 2) { warray[w_index++] = wr; warray[w_index++] = wi; double tempr = wr; wr = tempr * wfr - wi * wfi; wi = tempr * wfi + wi * wfr; } } // PRECOMPUTATION of wwr1, wwi1 for factor 4 Decomposition (3 * complex // operators and 8 +/- complex operators) { w_index = 0; int w_index2 = warray.length >> 1; for (int i = 0, nstep = 2; i < (imax - 1); i++) { int jmax = nstep; nstep *= 2; int ii = w_index + jmax; for (int j = 0; j < jmax; j += 2) { double wr = warray[w_index++]; double wi = warray[w_index++]; double wr1 = warray[ii++]; double wi1 = warray[ii++]; warray[w_index2++] = wr * wr1 - wi * wi1; warray[w_index2++] = wr * wi1 + wi * wr1; } } } return warray; } private final static void calc(int fftFrameSize, double[] data, int sign, double[] w) { final int fftFrameSize2 = fftFrameSize << 1; int nstep = 2; if (nstep >= fftFrameSize2) return; int i = nstep - 2; if (sign == -1) calcF4F(fftFrameSize, data, i, nstep, w); else calcF4I(fftFrameSize, data, i, nstep, w); } private final static void calcF2E(int fftFrameSize, double[] data, int i, int nstep, double[] w) { int jmax = nstep; for (int n = 0; n < jmax; n += 2) { double wr = w[i++]; double wi = w[i++]; int m = n + jmax; double datam_r = data[m]; double datam_i = data[m + 1]; double datan_r = data[n]; double datan_i = data[n + 1]; double tempr = datam_r * wr - datam_i * wi; double tempi = datam_r * wi + datam_i * wr; data[m] = datan_r - tempr; data[m + 1] = datan_i - tempi; data[n] = datan_r + tempr; data[n + 1] = datan_i + tempi; } return; } // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/- // complex operators private final static void calcF4F(int fftFrameSize, double[] data, int i, int nstep, double[] w) { final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize; // Factor-4 Decomposition int w_len = w.length >> 1; while (nstep < fftFrameSize2) { if (nstep << 2 == fftFrameSize2) { // Goto Factor-4 Final Decomposition // calcF4E(data, i, nstep, -1, w); calcF4FE(fftFrameSize, data, i, nstep, w); return; } int jmax = nstep; int nnstep = nstep << 1; if (nnstep == fftFrameSize2) { // Factor-4 Decomposition not possible calcF2E(fftFrameSize, data, i, nstep, w); return; } nstep <<= 2; int ii = i + jmax; int iii = i + w_len; { i += 2; ii += 2; iii += 2; for (int n = 0; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r; double tempi = datam1_i; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r; double n2w1i = datan2_i; double m2ww1r = datam2_r; double m2ww1i = datam2_i; tempr = m2ww1r - n2w1r; tempi = m2ww1i - n2w1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } for (int j = 2; j < jmax; j += 2) { double wr = w[i++]; double wi = w[i++]; double wr1 = w[ii++]; double wi1 = w[ii++]; double wwr1 = w[iii++]; double wwi1 = w[iii++]; // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be // precomputed!!! // double wwi1 = wr * wi1 + wi * wr1; for (int n = j; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r * wr - datam1_i * wi; double tempi = datam1_r * wi + datam1_i * wr; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r * wr1 - datan2_i * wi1; double n2w1i = datan2_r * wi1 + datan2_i * wr1; double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1; double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1; tempr = m2ww1r - n2w1r; tempi = m2ww1i - n2w1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } i += jmax << 1; } calcF2E(fftFrameSize, data, i, nstep, w); } // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/- // complex operators private final static void calcF4I(int fftFrameSize, double[] data, int i, int nstep, double[] w) { final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize; // Factor-4 Decomposition int w_len = w.length >> 1; while (nstep < fftFrameSize2) { if (nstep << 2 == fftFrameSize2) { // Goto Factor-4 Final Decomposition // calcF4E(data, i, nstep, 1, w); calcF4IE(fftFrameSize, data, i, nstep, w); return; } int jmax = nstep; int nnstep = nstep << 1; if (nnstep == fftFrameSize2) { // Factor-4 Decomposition not possible calcF2E(fftFrameSize, data, i, nstep, w); return; } nstep <<= 2; int ii = i + jmax; int iii = i + w_len; { i += 2; ii += 2; iii += 2; for (int n = 0; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r; double tempi = datam1_i; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r; double n2w1i = datan2_i; double m2ww1r = datam2_r; double m2ww1i = datam2_i; tempr = n2w1r - m2ww1r; tempi = n2w1i - m2ww1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } for (int j = 2; j < jmax; j += 2) { double wr = w[i++]; double wi = w[i++]; double wr1 = w[ii++]; double wi1 = w[ii++]; double wwr1 = w[iii++]; double wwi1 = w[iii++]; // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be // precomputed!!! // double wwi1 = wr * wi1 + wi * wr1; for (int n = j; n < fftFrameSize2; n += nstep) { int m = n + jmax; double datam1_r = data[m]; double datam1_i = data[m + 1]; double datan1_r = data[n]; double datan1_i = data[n + 1]; n += nnstep; m += nnstep; double datam2_r = data[m]; double datam2_i = data[m + 1]; double datan2_r = data[n]; double datan2_i = data[n + 1]; double tempr = datam1_r * wr - datam1_i * wi; double tempi = datam1_r * wi + datam1_i * wr; datam1_r = datan1_r - tempr; datam1_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; double n2w1r = datan2_r * wr1 - datan2_i * wi1; double n2w1i = datan2_r * wi1 + datan2_i * wr1; double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1; double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1; tempr = n2w1r - m2ww1r; tempi = n2w1i - m2ww1i; datam2_r = datam1_r + tempi; datam2_i = datam1_i - tempr; datam1_r = datam1_r - tempi; datam1_i = datam1_i + tempr; tempr = n2w1r + m2ww1r; tempi = n2w1i + m2ww1i; datan2_r = datan1_r - tempr; datan2_i = datan1_i - tempi; datan1_r = datan1_r + tempr; datan1_i = datan1_i + tempi; data[m] = datam2_r; data[m + 1] = datam2_i; data[n] = datan2_r; data[n + 1] = datan2_i; n -= nnstep; m -= nnstep; data[m] = datam1_r; data[m + 1] = datam1_i; data[n] = datan1_r; data[n + 1] = datan1_i; } } i += jmax << 1; } calcF2E(fftFrameSize, data, i, nstep, w); } // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-