1 /*
2  * Copyright (c) 2007, 2013, Oracle and/or its affiliates. All rights reserved.
3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4  *
5  * This code is free software; you can redistribute it and/or modify it
6  * under the terms of the GNU General Public License version 2 only, as
7  * published by the Free Software Foundation.  Oracle designates this
8  * particular file as subject to the "Classpath" exception as provided
9  * by Oracle in the LICENSE file that accompanied this code.
10  *
11  * This code is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14  * version 2 for more details (a copy is included in the LICENSE file that
15  * accompanied this code).
16  *
17  * You should have received a copy of the GNU General Public License version
18  * 2 along with this work; if not, write to the Free Software Foundation,
19  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20  *
21  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22  * or visit www.oracle.com if you need additional information or have any
23  * questions.
24  */
25 
26 package com.sun.media.sound;
27 
28 /**
29  * Fast Fourier Transformer.
30  *
31  * @author Karl Helgason
32  */
33 public final class FFT {
34 
35     private final double[] w;
36     private final int fftFrameSize;
37     private final int sign;
38     private final int[] bitm_array;
39     private final int fftFrameSize2;
40 
41     // Sign = -1 is FFT, 1 is IFFT (inverse FFT)
42     // Data = Interlaced double array to be transformed.
43     // The order is: real (sin), complex (cos)
44     // Framesize must be power of 2
FFT(int fftFrameSize, int sign)45     public FFT(int fftFrameSize, int sign) {
46         w = computeTwiddleFactors(fftFrameSize, sign);
47 
48         this.fftFrameSize = fftFrameSize;
49         this.sign = sign;
50         fftFrameSize2 = fftFrameSize << 1;
51 
52         // Pre-process Bit-Reversal
53         bitm_array = new int[fftFrameSize2];
54         for (int i = 2; i < fftFrameSize2; i += 2) {
55             int j;
56             int bitm;
57             for (bitm = 2, j = 0; bitm < fftFrameSize2; bitm <<= 1) {
58                 if ((i & bitm) != 0)
59                     j++;
60                 j <<= 1;
61             }
62             bitm_array[i] = j;
63         }
64 
65     }
66 
transform(double[] data)67     public void transform(double[] data) {
68         bitreversal(data);
69         calc(fftFrameSize, data, sign, w);
70     }
71 
computeTwiddleFactors(int fftFrameSize, int sign)72     private static double[] computeTwiddleFactors(int fftFrameSize,
73             int sign) {
74 
75         int imax = (int) (Math.log(fftFrameSize) / Math.log(2.));
76 
77         double[] warray = new double[(fftFrameSize - 1) * 4];
78         int w_index = 0;
79 
80         for (int i = 0,  nstep = 2; i < imax; i++) {
81             int jmax = nstep;
82             nstep <<= 1;
83 
84             double wr = 1.0;
85             double wi = 0.0;
86 
87             double arg = Math.PI / (jmax >> 1);
88             double wfr = Math.cos(arg);
89             double wfi = sign * Math.sin(arg);
90 
91             for (int j = 0; j < jmax; j += 2) {
92                 warray[w_index++] = wr;
93                 warray[w_index++] = wi;
94 
95                 double tempr = wr;
96                 wr = tempr * wfr - wi * wfi;
97                 wi = tempr * wfi + wi * wfr;
98             }
99         }
100 
101         // PRECOMPUTATION of wwr1, wwi1 for factor 4 Decomposition (3 * complex
102         // operators and 8 +/- complex operators)
103         {
104             w_index = 0;
105             int w_index2 = warray.length >> 1;
106             for (int i = 0,  nstep = 2; i < (imax - 1); i++) {
107                 int jmax = nstep;
108                 nstep *= 2;
109 
110                 int ii = w_index + jmax;
111                 for (int j = 0; j < jmax; j += 2) {
112                     double wr = warray[w_index++];
113                     double wi = warray[w_index++];
114                     double wr1 = warray[ii++];
115                     double wi1 = warray[ii++];
116                     warray[w_index2++] = wr * wr1 - wi * wi1;
117                     warray[w_index2++] = wr * wi1 + wi * wr1;
118                 }
119             }
120 
121         }
122 
123         return warray;
124     }
125 
calc(int fftFrameSize, double[] data, int sign, double[] w)126     private static void calc(int fftFrameSize, double[] data, int sign,
127             double[] w) {
128 
129         final int fftFrameSize2 = fftFrameSize << 1;
130 
131         int nstep = 2;
132 
133         if (nstep >= fftFrameSize2)
134             return;
135         int i = nstep - 2;
136         if (sign == -1)
137             calcF4F(fftFrameSize, data, i, nstep, w);
138         else
139             calcF4I(fftFrameSize, data, i, nstep, w);
140 
141     }
142 
calcF2E(int fftFrameSize, double[] data, int i, int nstep, double[] w)143     private static void calcF2E(int fftFrameSize, double[] data, int i,
144             int nstep, double[] w) {
145         int jmax = nstep;
146         for (int n = 0; n < jmax; n += 2) {
147             double wr = w[i++];
148             double wi = w[i++];
149             int m = n + jmax;
150             double datam_r = data[m];
151             double datam_i = data[m + 1];
152             double datan_r = data[n];
153             double datan_i = data[n + 1];
154             double tempr = datam_r * wr - datam_i * wi;
155             double tempi = datam_r * wi + datam_i * wr;
156             data[m] = datan_r - tempr;
157             data[m + 1] = datan_i - tempi;
158             data[n] = datan_r + tempr;
159             data[n + 1] = datan_i + tempi;
160         }
161         return;
162 
163     }
164 
165     // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
166     // complex operators
calcF4F(int fftFrameSize, double[] data, int i, int nstep, double[] w)167     private static void calcF4F(int fftFrameSize, double[] data, int i,
168             int nstep, double[] w) {
169         final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
170         // Factor-4 Decomposition
171 
172         int w_len = w.length >> 1;
173         while (nstep < fftFrameSize2) {
174 
175             if (nstep << 2 == fftFrameSize2) {
176                 // Goto Factor-4 Final Decomposition
177                 // calcF4E(data, i, nstep, -1, w);
178                 calcF4FE(fftFrameSize, data, i, nstep, w);
179                 return;
180             }
181             int jmax = nstep;
182             int nnstep = nstep << 1;
183             if (nnstep == fftFrameSize2) {
184                 // Factor-4 Decomposition not possible
185                 calcF2E(fftFrameSize, data, i, nstep, w);
186                 return;
187             }
188             nstep <<= 2;
189             int ii = i + jmax;
190             int iii = i + w_len;
191 
192             {
193                 i += 2;
194                 ii += 2;
195                 iii += 2;
196 
197                 for (int n = 0; n < fftFrameSize2; n += nstep) {
198                     int m = n + jmax;
199 
200                     double datam1_r = data[m];
201                     double datam1_i = data[m + 1];
202                     double datan1_r = data[n];
203                     double datan1_i = data[n + 1];
204 
205                     n += nnstep;
206                     m += nnstep;
207                     double datam2_r = data[m];
208                     double datam2_i = data[m + 1];
209                     double datan2_r = data[n];
210                     double datan2_i = data[n + 1];
211 
212                     double tempr = datam1_r;
213                     double tempi = datam1_i;
214 
215                     datam1_r = datan1_r - tempr;
216                     datam1_i = datan1_i - tempi;
217                     datan1_r = datan1_r + tempr;
218                     datan1_i = datan1_i + tempi;
219 
220                     double n2w1r = datan2_r;
221                     double n2w1i = datan2_i;
222                     double m2ww1r = datam2_r;
223                     double m2ww1i = datam2_i;
224 
225                     tempr = m2ww1r - n2w1r;
226                     tempi = m2ww1i - n2w1i;
227 
228                     datam2_r = datam1_r + tempi;
229                     datam2_i = datam1_i - tempr;
230                     datam1_r = datam1_r - tempi;
231                     datam1_i = datam1_i + tempr;
232 
233                     tempr = n2w1r + m2ww1r;
234                     tempi = n2w1i + m2ww1i;
235 
236                     datan2_r = datan1_r - tempr;
237                     datan2_i = datan1_i - tempi;
238                     datan1_r = datan1_r + tempr;
239                     datan1_i = datan1_i + tempi;
240 
241                     data[m] = datam2_r;
242                     data[m + 1] = datam2_i;
243                     data[n] = datan2_r;
244                     data[n + 1] = datan2_i;
245 
246                     n -= nnstep;
247                     m -= nnstep;
248                     data[m] = datam1_r;
249                     data[m + 1] = datam1_i;
250                     data[n] = datan1_r;
251                     data[n + 1] = datan1_i;
252 
253                 }
254             }
255 
256             for (int j = 2; j < jmax; j += 2) {
257                 double wr = w[i++];
258                 double wi = w[i++];
259                 double wr1 = w[ii++];
260                 double wi1 = w[ii++];
261                 double wwr1 = w[iii++];
262                 double wwi1 = w[iii++];
263                 // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
264                 // precomputed!!!
265                 // double wwi1 = wr * wi1 + wi * wr1;
266 
267                 for (int n = j; n < fftFrameSize2; n += nstep) {
268                     int m = n + jmax;
269 
270                     double datam1_r = data[m];
271                     double datam1_i = data[m + 1];
272                     double datan1_r = data[n];
273                     double datan1_i = data[n + 1];
274 
275                     n += nnstep;
276                     m += nnstep;
277                     double datam2_r = data[m];
278                     double datam2_i = data[m + 1];
279                     double datan2_r = data[n];
280                     double datan2_i = data[n + 1];
281 
282                     double tempr = datam1_r * wr - datam1_i * wi;
283                     double tempi = datam1_r * wi + datam1_i * wr;
284 
285                     datam1_r = datan1_r - tempr;
286                     datam1_i = datan1_i - tempi;
287                     datan1_r = datan1_r + tempr;
288                     datan1_i = datan1_i + tempi;
289 
290                     double n2w1r = datan2_r * wr1 - datan2_i * wi1;
291                     double n2w1i = datan2_r * wi1 + datan2_i * wr1;
292                     double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
293                     double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;
294 
295                     tempr = m2ww1r - n2w1r;
296                     tempi = m2ww1i - n2w1i;
297 
298                     datam2_r = datam1_r + tempi;
299                     datam2_i = datam1_i - tempr;
300                     datam1_r = datam1_r - tempi;
301                     datam1_i = datam1_i + tempr;
302 
303                     tempr = n2w1r + m2ww1r;
304                     tempi = n2w1i + m2ww1i;
305 
306                     datan2_r = datan1_r - tempr;
307                     datan2_i = datan1_i - tempi;
308                     datan1_r = datan1_r + tempr;
309                     datan1_i = datan1_i + tempi;
310 
311                     data[m] = datam2_r;
312                     data[m + 1] = datam2_i;
313                     data[n] = datan2_r;
314                     data[n + 1] = datan2_i;
315 
316                     n -= nnstep;
317                     m -= nnstep;
318                     data[m] = datam1_r;
319                     data[m + 1] = datam1_i;
320                     data[n] = datan1_r;
321                     data[n + 1] = datan1_i;
322                 }
323             }
324 
325             i += jmax << 1;
326 
327         }
328 
329         calcF2E(fftFrameSize, data, i, nstep, w);
330 
331     }
332 
333     // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
334     // complex operators
calcF4I(int fftFrameSize, double[] data, int i, int nstep, double[] w)335     private static void calcF4I(int fftFrameSize, double[] data, int i,
336             int nstep, double[] w) {
337         final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
338         // Factor-4 Decomposition
339 
340         int w_len = w.length >> 1;
341         while (nstep < fftFrameSize2) {
342 
343             if (nstep << 2 == fftFrameSize2) {
344                 // Goto Factor-4 Final Decomposition
345                 // calcF4E(data, i, nstep, 1, w);
346                 calcF4IE(fftFrameSize, data, i, nstep, w);
347                 return;
348             }
349             int jmax = nstep;
350             int nnstep = nstep << 1;
351             if (nnstep == fftFrameSize2) {
352                 // Factor-4 Decomposition not possible
353                 calcF2E(fftFrameSize, data, i, nstep, w);
354                 return;
355             }
356             nstep <<= 2;
357             int ii = i + jmax;
358             int iii = i + w_len;
359             {
360                 i += 2;
361                 ii += 2;
362                 iii += 2;
363 
364                 for (int n = 0; n < fftFrameSize2; n += nstep) {
365                     int m = n + jmax;
366 
367                     double datam1_r = data[m];
368                     double datam1_i = data[m + 1];
369                     double datan1_r = data[n];
370                     double datan1_i = data[n + 1];
371 
372                     n += nnstep;
373                     m += nnstep;
374                     double datam2_r = data[m];
375                     double datam2_i = data[m + 1];
376                     double datan2_r = data[n];
377                     double datan2_i = data[n + 1];
378 
379                     double tempr = datam1_r;
380                     double tempi = datam1_i;
381 
382                     datam1_r = datan1_r - tempr;
383                     datam1_i = datan1_i - tempi;
384                     datan1_r = datan1_r + tempr;
385                     datan1_i = datan1_i + tempi;
386 
387                     double n2w1r = datan2_r;
388                     double n2w1i = datan2_i;
389                     double m2ww1r = datam2_r;
390                     double m2ww1i = datam2_i;
391 
392                     tempr = n2w1r - m2ww1r;
393                     tempi = n2w1i - m2ww1i;
394 
395                     datam2_r = datam1_r + tempi;
396                     datam2_i = datam1_i - tempr;
397                     datam1_r = datam1_r - tempi;
398                     datam1_i = datam1_i + tempr;
399 
400                     tempr = n2w1r + m2ww1r;
401                     tempi = n2w1i + m2ww1i;
402 
403                     datan2_r = datan1_r - tempr;
404                     datan2_i = datan1_i - tempi;
405                     datan1_r = datan1_r + tempr;
406                     datan1_i = datan1_i + tempi;
407 
408                     data[m] = datam2_r;
409                     data[m + 1] = datam2_i;
410                     data[n] = datan2_r;
411                     data[n + 1] = datan2_i;
412 
413                     n -= nnstep;
414                     m -= nnstep;
415                     data[m] = datam1_r;
416                     data[m + 1] = datam1_i;
417                     data[n] = datan1_r;
418                     data[n + 1] = datan1_i;
419 
420                 }
421 
422             }
423             for (int j = 2; j < jmax; j += 2) {
424                 double wr = w[i++];
425                 double wi = w[i++];
426                 double wr1 = w[ii++];
427                 double wi1 = w[ii++];
428                 double wwr1 = w[iii++];
429                 double wwi1 = w[iii++];
430                 // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
431                 // precomputed!!!
432                 // double wwi1 = wr * wi1 + wi * wr1;
433 
434                 for (int n = j; n < fftFrameSize2; n += nstep) {
435                     int m = n + jmax;
436 
437                     double datam1_r = data[m];
438                     double datam1_i = data[m + 1];
439                     double datan1_r = data[n];
440                     double datan1_i = data[n + 1];
441 
442                     n += nnstep;
443                     m += nnstep;
444                     double datam2_r = data[m];
445                     double datam2_i = data[m + 1];
446                     double datan2_r = data[n];
447                     double datan2_i = data[n + 1];
448 
449                     double tempr = datam1_r * wr - datam1_i * wi;
450                     double tempi = datam1_r * wi + datam1_i * wr;
451 
452                     datam1_r = datan1_r - tempr;
453                     datam1_i = datan1_i - tempi;
454                     datan1_r = datan1_r + tempr;
455                     datan1_i = datan1_i + tempi;
456 
457                     double n2w1r = datan2_r * wr1 - datan2_i * wi1;
458                     double n2w1i = datan2_r * wi1 + datan2_i * wr1;
459                     double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
460                     double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;
461 
462                     tempr = n2w1r - m2ww1r;
463                     tempi = n2w1i - m2ww1i;
464 
465                     datam2_r = datam1_r + tempi;
466                     datam2_i = datam1_i - tempr;
467                     datam1_r = datam1_r - tempi;
468                     datam1_i = datam1_i + tempr;
469 
470                     tempr = n2w1r + m2ww1r;
471                     tempi = n2w1i + m2ww1i;
472 
473                     datan2_r = datan1_r - tempr;
474                     datan2_i = datan1_i - tempi;
475                     datan1_r = datan1_r + tempr;
476                     datan1_i = datan1_i + tempi;
477 
478                     data[m] = datam2_r;
479                     data[m + 1] = datam2_i;
480                     data[n] = datan2_r;
481                     data[n + 1] = datan2_i;
482 
483                     n -= nnstep;
484                     m -= nnstep;
485                     data[m] = datam1_r;
486                     data[m + 1] = datam1_i;
487                     data[n] = datan1_r;
488                     data[n + 1] = datan1_i;
489 
490                 }
491             }
492 
493             i += jmax << 1;
494 
495         }
496 
497         calcF2E(fftFrameSize, data, i, nstep, w);
498 
499     }
500 
501     // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
502     // complex operators
calcF4FE(int fftFrameSize, double[] data, int i, int nstep, double[] w)503     private static void calcF4FE(int fftFrameSize, double[] data, int i,
504             int nstep, double[] w) {
505         final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
506         // Factor-4 Decomposition
507 
508         int w_len = w.length >> 1;
509         while (nstep < fftFrameSize2) {
510 
511             int jmax = nstep;
512             int nnstep = nstep << 1;
513             if (nnstep == fftFrameSize2) {
514                 // Factor-4 Decomposition not possible
515                 calcF2E(fftFrameSize, data, i, nstep, w);
516                 return;
517             }
518             nstep <<= 2;
519             int ii = i + jmax;
520             int iii = i + w_len;
521             for (int n = 0; n < jmax; n += 2) {
522                 double wr = w[i++];
523                 double wi = w[i++];
524                 double wr1 = w[ii++];
525                 double wi1 = w[ii++];
526                 double wwr1 = w[iii++];
527                 double wwi1 = w[iii++];
528                 // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
529                 // precomputed!!!
530                 // double wwi1 = wr * wi1 + wi * wr1;
531 
532                 int m = n + jmax;
533 
534                 double datam1_r = data[m];
535                 double datam1_i = data[m + 1];
536                 double datan1_r = data[n];
537                 double datan1_i = data[n + 1];
538 
539                 n += nnstep;
540                 m += nnstep;
541                 double datam2_r = data[m];
542                 double datam2_i = data[m + 1];
543                 double datan2_r = data[n];
544                 double datan2_i = data[n + 1];
545 
546                 double tempr = datam1_r * wr - datam1_i * wi;
547                 double tempi = datam1_r * wi + datam1_i * wr;
548 
549                 datam1_r = datan1_r - tempr;
550                 datam1_i = datan1_i - tempi;
551                 datan1_r = datan1_r + tempr;
552                 datan1_i = datan1_i + tempi;
553 
554                 double n2w1r = datan2_r * wr1 - datan2_i * wi1;
555                 double n2w1i = datan2_r * wi1 + datan2_i * wr1;
556                 double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
557                 double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;
558 
559                 tempr = m2ww1r - n2w1r;
560                 tempi = m2ww1i - n2w1i;
561 
562                 datam2_r = datam1_r + tempi;
563                 datam2_i = datam1_i - tempr;
564                 datam1_r = datam1_r - tempi;
565                 datam1_i = datam1_i + tempr;
566 
567                 tempr = n2w1r + m2ww1r;
568                 tempi = n2w1i + m2ww1i;
569 
570                 datan2_r = datan1_r - tempr;
571                 datan2_i = datan1_i - tempi;
572                 datan1_r = datan1_r + tempr;
573                 datan1_i = datan1_i + tempi;
574 
575                 data[m] = datam2_r;
576                 data[m + 1] = datam2_i;
577                 data[n] = datan2_r;
578                 data[n + 1] = datan2_i;
579 
580                 n -= nnstep;
581                 m -= nnstep;
582                 data[m] = datam1_r;
583                 data[m + 1] = datam1_i;
584                 data[n] = datan1_r;
585                 data[n + 1] = datan1_i;
586 
587             }
588 
589             i += jmax << 1;
590 
591         }
592     }
593 
594     // Perform Factor-4 Decomposition with 3 * complex operators and 8 +/-
595     // complex operators
calcF4IE(int fftFrameSize, double[] data, int i, int nstep, double[] w)596     private static void calcF4IE(int fftFrameSize, double[] data, int i,
597             int nstep, double[] w) {
598         final int fftFrameSize2 = fftFrameSize << 1; // 2*fftFrameSize;
599         // Factor-4 Decomposition
600 
601         int w_len = w.length >> 1;
602         while (nstep < fftFrameSize2) {
603 
604             int jmax = nstep;
605             int nnstep = nstep << 1;
606             if (nnstep == fftFrameSize2) {
607                 // Factor-4 Decomposition not possible
608                 calcF2E(fftFrameSize, data, i, nstep, w);
609                 return;
610             }
611             nstep <<= 2;
612             int ii = i + jmax;
613             int iii = i + w_len;
614             for (int n = 0; n < jmax; n += 2) {
615                 double wr = w[i++];
616                 double wi = w[i++];
617                 double wr1 = w[ii++];
618                 double wi1 = w[ii++];
619                 double wwr1 = w[iii++];
620                 double wwi1 = w[iii++];
621                 // double wwr1 = wr * wr1 - wi * wi1; // these numbers can be
622                 // precomputed!!!
623                 // double wwi1 = wr * wi1 + wi * wr1;
624 
625                 int m = n + jmax;
626 
627                 double datam1_r = data[m];
628                 double datam1_i = data[m + 1];
629                 double datan1_r = data[n];
630                 double datan1_i = data[n + 1];
631 
632                 n += nnstep;
633                 m += nnstep;
634                 double datam2_r = data[m];
635                 double datam2_i = data[m + 1];
636                 double datan2_r = data[n];
637                 double datan2_i = data[n + 1];
638 
639                 double tempr = datam1_r * wr - datam1_i * wi;
640                 double tempi = datam1_r * wi + datam1_i * wr;
641 
642                 datam1_r = datan1_r - tempr;
643                 datam1_i = datan1_i - tempi;
644                 datan1_r = datan1_r + tempr;
645                 datan1_i = datan1_i + tempi;
646 
647                 double n2w1r = datan2_r * wr1 - datan2_i * wi1;
648                 double n2w1i = datan2_r * wi1 + datan2_i * wr1;
649                 double m2ww1r = datam2_r * wwr1 - datam2_i * wwi1;
650                 double m2ww1i = datam2_r * wwi1 + datam2_i * wwr1;
651 
652                 tempr = n2w1r - m2ww1r;
653                 tempi = n2w1i - m2ww1i;
654 
655                 datam2_r = datam1_r + tempi;
656                 datam2_i = datam1_i - tempr;
657                 datam1_r = datam1_r - tempi;
658                 datam1_i = datam1_i + tempr;
659 
660                 tempr = n2w1r + m2ww1r;
661                 tempi = n2w1i + m2ww1i;
662 
663                 datan2_r = datan1_r - tempr;
664                 datan2_i = datan1_i - tempi;
665                 datan1_r = datan1_r + tempr;
666                 datan1_i = datan1_i + tempi;
667 
668                 data[m] = datam2_r;
669                 data[m + 1] = datam2_i;
670                 data[n] = datan2_r;
671                 data[n + 1] = datan2_i;
672 
673                 n -= nnstep;
674                 m -= nnstep;
675                 data[m] = datam1_r;
676                 data[m + 1] = datam1_i;
677                 data[n] = datan1_r;
678                 data[n + 1] = datan1_i;
679 
680             }
681 
682             i += jmax << 1;
683 
684         }
685     }
686 
bitreversal(double[] data)687     private void bitreversal(double[] data) {
688         if (fftFrameSize < 4)
689             return;
690 
691         int inverse = fftFrameSize2 - 2;
692         for (int i = 0; i < fftFrameSize; i += 4) {
693             int j = bitm_array[i];
694 
695             // Performing Bit-Reversal, even v.s. even, O(2N)
696             if (i < j) {
697 
698                 int n = i;
699                 int m = j;
700 
701                 // COMPLEX: SWAP(data[n], data[m])
702                 // Real Part
703                 double tempr = data[n];
704                 data[n] = data[m];
705                 data[m] = tempr;
706                 // Imagery Part
707                 n++;
708                 m++;
709                 double tempi = data[n];
710                 data[n] = data[m];
711                 data[m] = tempi;
712 
713                 n = inverse - i;
714                 m = inverse - j;
715 
716                 // COMPLEX: SWAP(data[n], data[m])
717                 // Real Part
718                 tempr = data[n];
719                 data[n] = data[m];
720                 data[m] = tempr;
721                 // Imagery Part
722                 n++;
723                 m++;
724                 tempi = data[n];
725                 data[n] = data[m];
726                 data[m] = tempi;
727             }
728 
729             // Performing Bit-Reversal, odd v.s. even, O(N)
730 
731             int m = j + fftFrameSize; // bitm_array[i+2];
732             // COMPLEX: SWAP(data[n], data[m])
733             // Real Part
734             int n = i + 2;
735             double tempr = data[n];
736             data[n] = data[m];
737             data[m] = tempr;
738             // Imagery Part
739             n++;
740             m++;
741             double tempi = data[n];
742             data[n] = data[m];
743             data[m] = tempi;
744         }
745     }
746 }
747