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