2 //Title: 1-d mixed radix FFT.
4 //Copyright: Copyright (c) 1998
6 //Company: University of Wisconsin-Milwaukee.
8 // The number of DFT is factorized.
10 // Some short FFTs, such as length 2, 3, 4, 5, 8, 10, are used
11 // to improve the speed.
13 // Prime factors are processed using DFT. In the future, we can
15 // Note: there is no limit how large the prime factor can be,
16 // because for a set of data of an image, the length can be
17 // random, ie. an image can have size 263 x 300, where 263 is
18 // a large prime factor.
20 // A permute() function is used to make sure FFT can be calculated
23 // A triddle() function is used to perform the FFT.
25 // This program is for FFT of complex data, if the input is real,
26 // the program can be further improved. Because I want to use the
27 // same program to do IFFT, whose input is often complex, so I
28 // still use this program.
30 // To save the memory and improve the speed, float data are used
31 // instead of double, but I do have a double version transforms.fft.
33 // Factorize() is done in constructor, transforms.fft() is needed to be
34 // called to do FFT, this is good for use in fft2d, then
35 // factorize() is not needed for each row/column of data, since
36 // each row/column of a matrix has the same length.
41 // Maximum numbers of factors allowed.
42 //private int MaxFactorsNumber = 30;
43 public int MaxFactorsNumber;
45 // cos2to3PI = cos(2*pi/3), using for 3 point FFT.
46 // cos(2*PI/3) is not -1.5
47 public double cos2to3PI;
48 // sin2to3PI = sin(2*pi/3), using for 3 point FFT.
49 public double sin2to3PI;
51 // TwotoFivePI = 2*pi/5.
52 // c51, c52, c53, c54, c55 are used in fft5().
53 // c51 =(cos(TwotoFivePI)+cos(2*TwotoFivePI))/2-1.
55 // c52 =(cos(TwotoFivePI)-cos(2*TwotoFivePI))/2.
57 // c53 = -sin(TwotoFivePI).
59 // c54 =-(sin(TwotoFivePI)+sin(2*TwotoFivePI)).
61 // c55 =(sin(TwotoFivePI)-sin(2*TwotoFivePI)).
64 // OnetoSqrt2 = 1/sqrt(2), used in fft8().
65 public double OnetoSqrt2;
69 int N; // length of N point FFT.
70 int NumofFactors; // Number of factors of N.
71 int maxFactor; // Maximum factor of N.
73 int factors[]; // Factors of N processed in the current stage.
74 int sofar[]; // Finished factors before the current stage.
75 int remain[]; // Finished factors after the current stage.
77 double inputRe[], inputIm[]; // Input of FFT.
78 double temRe[], temIm[]; // Intermediate result of FFT.
79 double outputRe[], outputIm[]; // Output of FFT.
80 boolean factorsWerePrinted;
82 // Constructor: FFT of Complex data.
85 MaxFactorsNumber = 37;
87 sin2to3PI = 8.6602540378444E-01f;
89 c52 = 5.5901699437495E-01f;
90 c53 = -9.5105651629515E-01f;
91 c54 = -1.5388417685876E+00f;
92 c55 = 3.6327126400268E-01f;
93 OnetoSqrt2 = 7.0710678118655E-01f;
96 factorsWerePrinted = false;
97 outputRe = global new double[N];
98 outputIm = global new double[N];
103 // Allocate memory for intermediate result of FFT.
104 temRe = global new double[maxFactor]; //Check usage of this
105 temIm = global new double[maxFactor];
109 public void fft(double inputRe[], double inputIm[]) {
110 // First make sure inputRe & inputIm are of the same length.
111 if (inputRe.length != N || inputIm.length != N) {
112 System.printString("Error: the length of real part & imaginary part " +
113 "of the input to 1-d FFT are different");
116 this.inputRe = inputRe;
117 this.inputIm = inputIm;
120 //System.printString("ready to twiddle");
122 for (int factorIndex = 0; factorIndex < NumofFactors; factorIndex++)
123 twiddle(factorIndex);
124 //System.printString("ready to copy");
126 // Copy the output[] data to input[], so the output can be
127 // returned in the input array.
128 for (int i = 0; i < N; i++) {
129 inputRe[i] = outputRe[i];
130 inputIm[i] = outputIm[i];
136 public void printFactors() {
137 if (factorsWerePrinted) return;
138 factorsWerePrinted = true;
139 System.printString("factors.length = " + factors.length + "\n");
140 for (int i = 0; i < factors.length; i++)
141 System.printString("factors[i] = " + factors[i] + "\n");
144 public void factorize() {
145 int radices[] = global new int[6];
152 int temFactors[] = global new int[MaxFactorsNumber];
154 // 1 - point FFT, no need to factorize N.
160 // N - point FFT, N is needed to be factorized.
162 int index = 0; // index of temFactors.
163 int i = radices.length - 1;
165 while ((n > 1) && (i >= 0)) {
166 if ((n % radices[i]) == 0) {
168 temFactors[index++] = radices[i];
173 // Substitute 2x8 with 4x4.
174 // index>0, in the case only one prime factor, such as N=263.
175 if ((index > 0) && (temFactors[index - 1] == 2)) {
177 for (i = index - 2; (i >= 0) && (test == 0); i--) {
178 if (temFactors[i] == 8) {
179 temFactors[index - 1] = temFactors[i] = 4;
180 // break out of for loop, because only one '2' will exist in
181 // temFactors, so only one substitutation is needed.
189 for (int k = 2; k < Math.sqrt(n) + 1; k++)
190 while ((n % k) == 0) {
192 temFactors[index++] = k;
195 temFactors[index++] = n;
198 NumofFactors = index;
199 //if(temFactors[NumofFactors-1] > 10)
204 // Inverse temFactors and store factors into factors[].
205 factors = global new int[NumofFactors];
206 for (i = 0; i < NumofFactors; i++) {
207 factors[i] = temFactors[NumofFactors - i - 1];
210 // Calculate sofar[], remain[].
211 // sofar[] : finished factors before the current stage.
212 // factors[]: factors of N processed in the current stage.
213 // remain[] : finished factors after the current stage.
214 sofar = global new int[NumofFactors];
215 remain = global new int[NumofFactors];
217 remain[0] = N / factors[0];
219 for (i = 1; i < NumofFactors; i++) {
220 sofar[i] = sofar[i - 1] * factors[i - 1];
221 remain[i] = remain[i - 1] / factors[i];
223 } // End of function factorize().
225 private void permute() {
226 int count[] = new int[MaxFactorsNumber];
230 for (int i = 0; i < N - 1; i++) {
231 outputRe[i] = inputRe[k];
232 outputIm[i] = inputIm[k];
235 count[0] = count[0] + 1;
236 while (count[j] >= factors[j]) {
238 k = k - (j == 0?N:remain[j - 1]) + remain[j + 1];
240 count[j] = count[j] + 1;
243 outputRe[N - 1] = inputRe[N - 1];
244 outputIm[N - 1] = inputIm[N - 1];
245 } // End of function permute().
248 private void twiddle(int factorIndex) {
250 int sofarRadix = sofar[factorIndex];
251 int radix = factors[factorIndex];
252 int remainRadix = remain[factorIndex];
254 double tem; // Temporary variable to do data exchange.
256 double W = 2 * (double) Math.PI / (sofarRadix * radix);
257 double cosW = (double) Math.cos(W);
258 double sinW = -(double) Math.sin(W);
260 double twiddleRe[] = new double[radix];
261 double twiddleIm[] = new double[radix];
262 double twRe = 1.0f, twIm = 0f;
264 //Initialize twiddle addBk.address variables.
265 int dataOffset = 0, groupOffset = 0, address = 0;
267 for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
268 //System.printString("datano="+dataNo);
269 if (sofarRadix > 1) {
274 for (int i = 2; i < radix; i++) {
277 twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
278 twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
280 tem = cosW * twRe - sinW * twIm;
281 twIm = sinW * twRe + cosW * twIm;
284 for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
285 //System.printString("groupNo="+groupNo);
286 if ((sofarRadix > 1) && (dataNo > 0)) {
287 temRe[0] = outputRe[address];
288 temIm[0] = outputIm[address];
291 address = address + sofarRadix;
292 temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
293 twiddleIm[blockIndex] * outputIm[address];
294 temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
295 twiddleIm[blockIndex] * outputRe[address];
297 } while (blockIndex < radix);
299 for (int i = 0; i < radix; i++) {
300 //System.printString("temRe.length="+temRe.length);
301 //System.printString("i = "+i);
302 temRe[i] = outputRe[address];
303 temIm[i] = outputIm[address];
304 address += sofarRadix;
306 //System.printString("radix="+radix);
309 tem = temRe[0] + temRe[1];
310 temRe[1] = temRe[0] - temRe[1];
312 tem = temIm[0] + temIm[1];
313 temIm[1] = temIm[0] - temIm[1];
317 double t1Re = temRe[1] + temRe[2];
318 double t1Im = temIm[1] + temIm[2];
319 temRe[0] = temRe[0] + t1Re;
320 temIm[0] = temIm[0] + t1Im;
322 double m1Re = cos2to3PI * t1Re;
323 double m1Im = cos2to3PI * t1Im;
324 double m2Re = sin2to3PI * (temIm[1] - temIm[2]);
325 double m2Im = sin2to3PI * (temRe[2] - temRe[1]);
326 double s1Re = temRe[0] + m1Re;
327 double s1Im = temIm[0] + m1Im;
329 temRe[1] = s1Re + m2Re;
330 temIm[1] = s1Im + m2Im;
331 temRe[2] = s1Re - m2Re;
332 temIm[2] = s1Im - m2Im;
350 address = groupOffset;
351 for (int i = 0; i < radix; i++) {
352 outputRe[address] = temRe[i];
353 outputIm[address] = temIm[i];
354 address += sofarRadix;
356 groupOffset += sofarRadix * radix;
357 address = groupOffset;
359 groupOffset = ++dataOffset;
360 address = groupOffset;
362 } // End of function twiddle().
365 // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
366 private void fft4(double dataRe[], double dataIm[]) {
367 double t1Re,t1Im, t2Re,t2Im;
368 double m2Re,m2Im, m3Re,m3Im;
370 t1Re = dataRe[0] + dataRe[2];
371 t1Im = dataIm[0] + dataIm[2];
372 t2Re = dataRe[1] + dataRe[3];
373 t2Im = dataIm[1] + dataIm[3];
375 m2Re = dataRe[0] - dataRe[2];
376 m2Im = dataIm[0] - dataIm[2];
377 m3Re = dataIm[1] - dataIm[3];
378 m3Im = dataRe[3] - dataRe[1];
380 dataRe[0] = t1Re + t2Re;
381 dataIm[0] = t1Im + t2Im;
382 dataRe[2] = t1Re - t2Re;
383 dataIm[2] = t1Im - t2Im;
384 dataRe[1] = m2Re + m3Re;
385 dataIm[1] = m2Im + m3Im;
386 dataRe[3] = m2Re - m3Re;
387 dataIm[3] = m2Im - m3Im;
388 } // End of function fft4().
391 // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
392 private void fft5(double dataRe[], double dataIm[]) {
393 double t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
394 double m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
395 double s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
397 t1Re = dataRe[1] + dataRe[4];
398 t1Im = dataIm[1] + dataIm[4];
399 t2Re = dataRe[2] + dataRe[3];
400 t2Im = dataIm[2] + dataIm[3];
401 t3Re = dataRe[1] - dataRe[4];
402 t3Im = dataIm[1] - dataIm[4];
403 t4Re = dataRe[3] - dataRe[2];
404 t4Im = dataIm[3] - dataIm[2];
408 dataRe[0] = dataRe[0] + t5Re;
409 dataIm[0] = dataIm[0] + t5Im;
413 m2Re = c52 * (t1Re - t2Re);
414 m2Im = c52 * (t1Im - t2Im);
415 m3Re = -c53 * (t3Im + t4Im);
416 m3Im = c53 * (t3Re + t4Re);
426 s1Re = dataRe[0] + m1Re;
427 s1Im = dataIm[0] + m1Im;
433 dataRe[1] = s2Re + s3Re;
434 dataIm[1] = s2Im + s3Im;
435 dataRe[2] = s4Re + s5Re;
436 dataIm[2] = s4Im + s5Im;
437 dataRe[3] = s4Re - s5Re;
438 dataIm[3] = s4Im - s5Im;
439 dataRe[4] = s2Re - s3Re;
440 dataIm[4] = s2Im - s3Im;
441 } // End of function fft5().
445 private void fft8() {
446 double data1Re[] = new double[4];
447 double data1Im[] = new double[4];
448 double data2Re[] = new double[4];
449 double data2Im[] = new double[4];
452 // To improve the speed, use direct assaignment instead for loop here.
453 data1Re[0] = temRe[0];
454 data2Re[0] = temRe[1];
455 data1Re[1] = temRe[2];
456 data2Re[1] = temRe[3];
457 data1Re[2] = temRe[4];
458 data2Re[2] = temRe[5];
459 data1Re[3] = temRe[6];
460 data2Re[3] = temRe[7];
462 data1Im[0] = temIm[0];
463 data2Im[0] = temIm[1];
464 data1Im[1] = temIm[2];
465 data2Im[1] = temIm[3];
466 data1Im[2] = temIm[4];
467 data2Im[2] = temIm[5];
468 data1Im[3] = temIm[6];
469 data2Im[3] = temIm[7];
471 fft4(data1Re, data1Im);
472 fft4(data2Re, data2Im);
474 tem = OnetoSqrt2 * (data2Re[1] + data2Im[1]);
475 data2Im[1] = OnetoSqrt2 * (data2Im[1] - data2Re[1]);
478 data2Im[2] = -data2Re[2];
480 tem = OnetoSqrt2 * (data2Im[3] - data2Re[3]);
481 data2Im[3] = -OnetoSqrt2 * (data2Re[3] + data2Im[3]);
484 temRe[0] = data1Re[0] + data2Re[0];
485 temRe[4] = data1Re[0] - data2Re[0];
486 temRe[1] = data1Re[1] + data2Re[1];
487 temRe[5] = data1Re[1] - data2Re[1];
488 temRe[2] = data1Re[2] + data2Re[2];
489 temRe[6] = data1Re[2] - data2Re[2];
490 temRe[3] = data1Re[3] + data2Re[3];
491 temRe[7] = data1Re[3] - data2Re[3];
493 temIm[0] = data1Im[0] + data2Im[0];
494 temIm[4] = data1Im[0] - data2Im[0];
495 temIm[1] = data1Im[1] + data2Im[1];
496 temIm[5] = data1Im[1] - data2Im[1];
497 temIm[2] = data1Im[2] + data2Im[2];
498 temIm[6] = data1Im[2] - data2Im[2];
499 temIm[3] = data1Im[3] + data2Im[3];
500 temIm[7] = data1Im[3] - data2Im[3];
501 } // End of function fft8().
505 private void fft10() {
506 double data1Re[] = new double[5];
507 double data1Im[] = new double[5];
508 double data2Re[] = new double[5];
509 double data2Im[] = new double[5];
511 // To improve the speed, use direct assaignment instead for loop here.
512 data1Re[0] = temRe[0];
513 data2Re[0] = temRe[5];
514 data1Re[1] = temRe[2];
515 data2Re[1] = temRe[7];
516 data1Re[2] = temRe[4];
517 data2Re[2] = temRe[9];
518 data1Re[3] = temRe[6];
519 data2Re[3] = temRe[1];
520 data1Re[4] = temRe[8];
521 data2Re[4] = temRe[3];
523 data1Im[0] = temIm[0];
524 data2Im[0] = temIm[5];
525 data1Im[1] = temIm[2];
526 data2Im[1] = temIm[7];
527 data1Im[2] = temIm[4];
528 data2Im[2] = temIm[9];
529 data1Im[3] = temIm[6];
530 data2Im[3] = temIm[1];
531 data1Im[4] = temIm[8];
532 data2Im[4] = temIm[3];
534 fft5(data1Re, data1Im);
535 fft5(data2Re, data2Im);
537 temRe[0] = data1Re[0] + data2Re[0];
538 temRe[5] = data1Re[0] - data2Re[0];
539 temRe[6] = data1Re[1] + data2Re[1];
540 temRe[1] = data1Re[1] - data2Re[1];
541 temRe[2] = data1Re[2] + data2Re[2];
542 temRe[7] = data1Re[2] - data2Re[2];
543 temRe[8] = data1Re[3] + data2Re[3];
544 temRe[3] = data1Re[3] - data2Re[3];
545 temRe[4] = data1Re[4] + data2Re[4];
546 temRe[9] = data1Re[4] - data2Re[4];
548 temIm[0] = data1Im[0] + data2Im[0];
549 temIm[5] = data1Im[0] - data2Im[0];
550 temIm[6] = data1Im[1] + data2Im[1];
551 temIm[1] = data1Im[1] - data2Im[1];
552 temIm[2] = data1Im[2] + data2Im[2];
553 temIm[7] = data1Im[2] - data2Im[2];
554 temIm[8] = data1Im[3] + data2Im[3];
555 temIm[3] = data1Im[3] - data2Im[3];
556 temIm[4] = data1Im[4] + data2Im[4];
557 temIm[9] = data1Im[4] - data2Im[4];
558 } // End of function fft10().
562 public double sqrt(double d) {
568 private void fftPrime(int radix) {
570 double W = 2 * (double) Math.PI / radix;
571 double cosW = (double) Math.cos(W);
572 double sinW = -(double) Math.sin(W);
573 double WRe[] = new double[radix];
574 double WIm[] = new double[radix];
581 for (int i = 2; i < radix; i++) {
582 WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
583 WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
586 // FFT of prime length data, using DFT, can be improved in the future.
587 double rere, reim, imre, imim;
589 int max = (radix + 1) / 2;
591 double tem1Re[] = new double[max];
592 double tem1Im[] = new double[max];
593 double tem2Re[] = new double[max];
594 double tem2Im[] = new double[max];
596 for (j = 1; j < max; j++) {
597 tem1Re[j] = temRe[j] + temRe[radix - j];
598 tem1Im[j] = temIm[j] - temIm[radix - j];
599 tem2Re[j] = temRe[j] - temRe[radix - j];
600 tem2Im[j] = temIm[j] + temIm[radix - j];
603 for (j = 1; j < max; j++) {
606 temRe[radix - j] = temRe[0];
607 temIm[radix - j] = temIm[0];
609 for (int i = 1; i < max; i++) {
610 rere = WRe[k] * tem1Re[i];
611 imim = WIm[k] * tem1Im[i];
612 reim = WRe[k] * tem2Im[i];
613 imre = WIm[k] * tem2Re[i];
615 temRe[radix - j] += rere + imim;
616 temIm[radix - j] += reim - imre;
617 temRe[j] += rere - imim;
618 temIm[j] += reim + imre;
625 for (j = 1; j < max; j++) {
626 temRe[0] = temRe[0] + tem1Re[j];
627 temIm[0] = temIm[0] + tem2Im[j];
629 } // End of function fftPrime().
632 } // End of class FFT2d