|
/Users/lyon/j4p/src/math/transforms/fft/FFT1d.java
|
1 package math.transforms.fft;
2
3 //Title: 1-d mixed radix FFT.
4 //Version:
5 //Copyright: Copyright (c) 1998
6 //Author: Dongyan Wang
7 //Company: University of Wisconsin-Milwaukee.
8 //Description:
9 // The number of DFT is factorized.
10 //
11 // Some short FFTs, such as length 2, 3, 4, 5, 8, 10, are used
12 // to improve the speed.
13 //
14 // Prime factors are processed using DFT. In the future, we can
15 // improve this part.
16 // Note: there is no limit how large the prime factor can be,
17 // because for a set of data of an image, the length can be
18 // random, ie. an image can have size 263 x 300, where 263 is
19 // a large prime factor.
20 //
21 // A permute() function is used to make sure FFT can be calculated
22 // in place.
23 //
24 // A triddle() function is used to perform the FFT.
25 //
26 // This program is for FFT of complex data, if the input is real,
27 // the program can be further improved. Because I want to use the
28 // same program to do IFFT, whose input is often complex, so I
29 // still use this program.
30 //
31 // To save the memory and improve the speed, float data are used
32 // instead of double, but I do have a double version transforms.fft.
33 //
34 // Factorize() is done in constructor, transforms.fft() is needed to be
35 // called to do FFT, this is good for use in fft2d, then
36 // factorize() is not needed for each row/column of data, since
37 // each row/column of a matrix has the same length.
38 //
39
40
41 public class FFT1d extends TempVars {
42 // Maximum numbers of factors allowed.
43 //private static final int MaxFactorsNumber = 30;
44 private static final int MaxFactorsNumber = 37;
45
46 // cos2to3PI = cos(2*pi/3), using for 3 point FFT.
47 // cos(2*PI/3) is not -1.5
48 private static final float cos2to3PI = -1.5000f;
49 // sin2to3PI = sin(2*pi/3), using for 3 point FFT.
50 private static final float sin2to3PI = 8.6602540378444E-01f;
51
52 // TwotoFivePI = 2*pi/5.
53 // c51, c52, c53, c54, c55 are used in fft5().
54 // c51 =(cos(TwotoFivePI)+cos(2*TwotoFivePI))/2-1.
55 private static final float c51 = -1.25f;
56 // c52 =(cos(TwotoFivePI)-cos(2*TwotoFivePI))/2.
57 private static final float c52 = 5.5901699437495E-01f;
58 // c53 = -sin(TwotoFivePI).
59 private static final float c53 = -9.5105651629515E-01f;
60 // c54 =-(sin(TwotoFivePI)+sin(2*TwotoFivePI)).
61 private static final float c54 = -1.5388417685876E+00f;
62 // c55 =(sin(TwotoFivePI)-sin(2*TwotoFivePI)).
63 private static final float c55 = 3.6327126400268E-01f;
64
65 // OnetoSqrt2 = 1/sqrt(2), used in fft8().
66 private static final float OnetoSqrt2 = 7.0710678118655E-01f;
67
68 private static int lastRadix = 0;
69
70 int N; // length of N point FFT.
71 int NumofFactors; // Number of factors of N.
72 static final int maxFactor = 20; // Maximum factor of N.
73
74 int factors[]; // Factors of N processed in the current stage.
75 int sofar[]; // Finished factors before the current stage.
76 int remain[]; // Finished factors after the current stage.
77
78 float inputRe[], inputIm[]; // Input of FFT.
79 float temRe[], temIm[]; // Intermediate result of FFT.
80 float outputRe[], outputIm[]; // Output of FFT.
81 static boolean factorsWerePrinted = false;
82
83 // Constructor: FFT of Complex data.
84 public FFT1d(int N) {
85 this.N = N;
86 outputRe = new float[N];
87 outputIm = new float[N];
88
89 factorize();
90 printFactors();
91
92 // Allocate memory for intermediate result of FFT.
93 temRe = new float[maxFactor];
94 temIm = new float[maxFactor];
95 }
96
97 public void fft(float inputRe[], float inputIm[]) {
98 // First make sure inputRe & inputIm are of the same length.
99 if (inputRe.length != N || inputIm.length != N) {
100 System.out.println("Error: the length of real part & imaginary part " +
101 "of the input to 1-d FFT are different");
102 return;
103 } else {
104 this.inputRe = inputRe;
105 this.inputIm = inputIm;
106
107 permute();
108 //System.out.println("ready to twiddle");
109
110 for (int factorIndex = 0; factorIndex < NumofFactors; factorIndex++)
111 twiddle(factorIndex);
112 //System.out.println("ready to copy");
113
114 // Copy the output[] data to input[], so the output can be
115 // returned in the input array.
116 for (int i = 0; i < N; i++) {
117 inputRe[i] = outputRe[i];
118 inputIm[i] = outputIm[i];
119 }
120
121 }
122 }
123
124 public void printFactors() {
125 if (factorsWerePrinted) return;
126 factorsWerePrinted = true;
127 for (int i = 0; i < factors.length; i++)
128 System.out.println("factors[i] = " + factors[i]);
129 }
130
131 private void factorize() {
132 int radices[] = {2, 3, 4, 5, 8, 10};
133 int temFactors[] = new int[MaxFactorsNumber];
134
135 // 1 - point FFT, no need to factorize N.
136 if (N == 1) {
137 temFactors[0] = 1;
138 NumofFactors = 1;
139 }
140
141 // N - point FFT, N is needed to be factorized.
142 int n = N;
143 int index = 0; // index of temFactors.
144 int i = radices.length - 1;
145
146 while ((n > 1) && (i >= 0)) {
147 if ((n % radices[i]) == 0) {
148 n /= radices[i];
149 temFactors[index++] = radices[i];
150 } else
151 i--;
152 }
153
154 // Substitute 2x8 with 4x4.
155 // index>0, in the case only one prime factor, such as N=263.
156 if ((index > 0) && (temFactors[index - 1] == 2))
157 for (i = index - 2; i >= 0; i--)
158 if (temFactors[i] == 8) {
159 temFactors[index - 1] = temFactors[i] = 4;
160 // break out of for loop, because only one '2' will exist in
161 // temFactors, so only one substitutation is needed.
162 break;
163 }
164
165 if (n > 1) {
166 for (int k = 2; k < Math.sqrt(n) + 1; k++)
167 while ((n % k) == 0) {
168 n /= k;
169 temFactors[index++] = k;
170 }
171 if (n > 1) {
172 temFactors[index++] = n;
173 }
174 }
175 NumofFactors = index;
176 //if(temFactors[NumofFactors-1] > 10)
177 // maxFactor = n;
178 //else
179 // maxFactor = 10;
180
181 // Inverse temFactors and store factors into factors[].
182 factors = new int[NumofFactors];
183 for (i = 0; i < NumofFactors; i++) {
184 factors[i] = temFactors[NumofFactors - i - 1];
185 }
186
187 // Calculate sofar[], remain[].
188 // sofar[] : finished factors before the current stage.
189 // factors[]: factors of N processed in the current stage.
190 // remain[] : finished factors after the current stage.
191 sofar = new int[NumofFactors];
192 remain = new int[NumofFactors];
193
194 remain[0] = N / factors[0];
195 sofar[0] = 1;
196 for (i = 1; i < NumofFactors; i++) {
197 sofar[i] = sofar[i - 1] * factors[i - 1];
198 remain[i] = remain[i - 1] / factors[i];
199 }
200 } // End of function factorize().
201
202 private void permute() {
203 int count[] = new int[MaxFactorsNumber];
204 int j;
205 int k = 0;
206
207 for (int i = 0; i < N - 1; i++) {
208 outputRe[i] = inputRe[k];
209 outputIm[i] = inputIm[k];
210 j = 0;
211 k = k + remain[j];
212 count[0] = count[0] + 1;
213 while (count[j] >= factors[j]) {
214 count[j] = 0;
215 k = k - (j == 0?N:remain[j - 1]) + remain[j + 1];
216 j++;
217 count[j] = count[j] + 1;
218 }
219 }
220 outputRe[N - 1] = inputRe[N - 1];
221 outputIm[N - 1] = inputIm[N - 1];
222 } // End of function permute().
223
224 private void twiddle(int factorIndex) {
225 // Get factor data.
226 int sofarRadix = sofar[factorIndex];
227 int radix = factors[factorIndex];
228 int remainRadix = remain[factorIndex];
229
230 float tem; // Temporary variable to do data exchange.
231
232 float W = 2 * (float) Math.PI / (sofarRadix * radix);
233 float cosW = (float) Math.cos(W);
234 float sinW = -(float) Math.sin(W);
235
236 float twiddleRe[] = new float[radix];
237 float twiddleIm[] = new float[radix];
238 float twRe = 1.0f, twIm = 0f;
239
240 //Initialize twiddle addBk.address variables.
241 int dataOffset = 0, groupOffset = 0, address = 0;
242
243 for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
244 //System.out.println("datano="+dataNo);
245 if (sofarRadix > 1) {
246 twiddleRe[0] = 1.0f;
247 twiddleIm[0] = 0.0f;
248 twiddleRe[1] = twRe;
249 twiddleIm[1] = twIm;
250 for (int i = 2; i < radix; i++) {
251
252
253 twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
254 twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
255 }
256 tem = cosW * twRe - sinW * twIm;
257 twIm = sinW * twRe + cosW * twIm;
258 twRe = tem;
259 }
260 for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
261 //System.out.println("groupNo="+groupNo);
262 if ((sofarRadix > 1) && (dataNo > 0)) {
263 temRe[0] = outputRe[address];
264 temIm[0] = outputIm[address];
265 int blockIndex = 1;
266 do {
267 address = address + sofarRadix;
268 temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
269 twiddleIm[blockIndex] * outputIm[address];
270 temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
271 twiddleIm[blockIndex] * outputRe[address];
272 blockIndex++;
273 } while (blockIndex < radix);
274 } else
275 for (int i = 0; i < radix; i++) {
276 //System.out.println("temRe.length="+temRe.length);
277 //System.out.println("i = "+i);
278 temRe[i] = outputRe[address];
279 temIm[i] = outputIm[address];
280 address += sofarRadix;
281 }
282 //System.out.println("radix="+radix);
283 switch (radix) {
284
285
286 case 2:
287 tem = temRe[0] + temRe[1];
288 temRe[1] = temRe[0] - temRe[1];
289 temRe[0] = tem;
290 tem = temIm[0] + temIm[1];
291 temIm[1] = temIm[0] - temIm[1];
292 temIm[0] = tem;
293 break;
294 case 3:
295 float t1Re = temRe[1] + temRe[2];
296 float t1Im = temIm[1] + temIm[2];
297 temRe[0] = temRe[0] + t1Re;
298 temIm[0] = temIm[0] + t1Im;
299
300 float m1Re = cos2to3PI * t1Re;
301 float m1Im = cos2to3PI * t1Im;
302 float m2Re = sin2to3PI * (temIm[1] - temIm[2]);
303 float m2Im = sin2to3PI * (temRe[2] - temRe[1]);
304 float s1Re = temRe[0] + m1Re;
305 float s1Im = temIm[0] + m1Im;
306
307 temRe[1] = s1Re + m2Re;
308 temIm[1] = s1Im + m2Im;
309 temRe[2] = s1Re - m2Re;
310 temIm[2] = s1Im - m2Im;
311 break;
312 case 4:
313 fft4(temRe, temIm);
314 break;
315 case 5:
316 fft5(temRe, temIm);
317 break;
318 case 8:
319 fft8();
320 break;
321 case 10:
322 fft10();
323 break;
324 default :
325 fftPrime(radix);
326 break;
327 }
328 address = groupOffset;
329 for (int i = 0; i < radix; i++) {
330 outputRe[address] = temRe[i];
331 outputIm[address] = temIm[i];
332 address += sofarRadix;
333 }
334 groupOffset += sofarRadix * radix;
335 address = groupOffset;
336 }
337 groupOffset = ++dataOffset;
338 address = groupOffset;
339 }
340 } // End of function twiddle().
341
342 // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
343 private void fft4(float dataRe[], float dataIm[]) {
344 float t1Re,t1Im, t2Re,t2Im;
345 float m2Re,m2Im, m3Re,m3Im;
346
347 t1Re = dataRe[0] + dataRe[2];
348 t1Im = dataIm[0] + dataIm[2];
349 t2Re = dataRe[1] + dataRe[3];
350 t2Im = dataIm[1] + dataIm[3];
351
352 m2Re = dataRe[0] - dataRe[2];
353 m2Im = dataIm[0] - dataIm[2];
354 m3Re = dataIm[1] - dataIm[3];
355 m3Im = dataRe[3] - dataRe[1];
356
357 dataRe[0] = t1Re + t2Re;
358 dataIm[0] = t1Im + t2Im;
359 dataRe[2] = t1Re - t2Re;
360 dataIm[2] = t1Im - t2Im;
361 dataRe[1] = m2Re + m3Re;
362 dataIm[1] = m2Im + m3Im;
363 dataRe[3] = m2Re - m3Re;
364 dataIm[3] = m2Im - m3Im;
365 } // End of function fft4().
366
367 // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
368 private void fft5(float dataRe[], float dataIm[]) {
369 float t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
370 float m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
371 float s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
372
373 t1Re = dataRe[1] + dataRe[4];
374 t1Im = dataIm[1] + dataIm[4];
375 t2Re = dataRe[2] + dataRe[3];
376 t2Im = dataIm[2] + dataIm[3];
377 t3Re = dataRe[1] - dataRe[4];
378 t3Im = dataIm[1] - dataIm[4];
379 t4Re = dataRe[3] - dataRe[2];
380 t4Im = dataIm[3] - dataIm[2];
381 t5Re = t1Re + t2Re;
382 t5Im = t1Im + t2Im;
383
384 dataRe[0] = dataRe[0] + t5Re;
385 dataIm[0] = dataIm[0] + t5Im;
386
387 m1Re = c51 * t5Re;
388 m1Im = c51 * t5Im;
389 m2Re = c52 * (t1Re - t2Re);
390 m2Im = c52 * (t1Im - t2Im);
391 m3Re = -c53 * (t3Im + t4Im);
392 m3Im = c53 * (t3Re + t4Re);
393 m4Re = -c54 * t4Im;
394 m4Im = c54 * t4Re;
395 m5Re = -c55 * t3Im;
396 m5Im = c55 * t3Re;
397
398 s3Re = m3Re - m4Re;
399 s3Im = m3Im - m4Im;
400 s5Re = m3Re + m5Re;
401 s5Im = m3Im + m5Im;
402 s1Re = dataRe[0] + m1Re;
403 s1Im = dataIm[0] + m1Im;
404 s2Re = s1Re + m2Re;
405 s2Im = s1Im + m2Im;
406 s4Re = s1Re - m2Re;
407 s4Im = s1Im - m2Im;
408
409 dataRe[1] = s2Re + s3Re;
410 dataIm[1] = s2Im + s3Im;
411 dataRe[2] = s4Re + s5Re;
412 dataIm[2] = s4Im + s5Im;
413 dataRe[3] = s4Re - s5Re;
414 dataIm[3] = s4Im - s5Im;
415 dataRe[4] = s2Re - s3Re;
416 dataIm[4] = s2Im - s3Im;
417 } // End of function fft5().
418
419 private void fft8() {
420 float data1Re[] = new float[4];
421 float data1Im[] = new float[4];
422 float data2Re[] = new float[4];
423 float data2Im[] = new float[4];
424 float tem;
425
426 // To improve the speed, use direct assaignment instead for loop here.
427 data1Re[0] = temRe[0];
428 data2Re[0] = temRe[1];
429 data1Re[1] = temRe[2];
430 data2Re[1] = temRe[3];
431 data1Re[2] = temRe[4];
432 data2Re[2] = temRe[5];
433 data1Re[3] = temRe[6];
434 data2Re[3] = temRe[7];
435
436 data1Im[0] = temIm[0];
437 data2Im[0] = temIm[1];
438 data1Im[1] = temIm[2];
439 data2Im[1] = temIm[3];
440 data1Im[2] = temIm[4];
441 data2Im[2] = temIm[5];
442 data1Im[3] = temIm[6];
443 data2Im[3] = temIm[7];
444
445 fft4(data1Re, data1Im);
446 fft4(data2Re, data2Im);
447
448 tem = OnetoSqrt2 * (data2Re[1] + data2Im[1]);
449 data2Im[1] = OnetoSqrt2 * (data2Im[1] - data2Re[1]);
450 data2Re[1] = tem;
451 tem = data2Im[2];
452 data2Im[2] = -data2Re[2];
453 data2Re[2] = tem;
454 tem = OnetoSqrt2 * (data2Im[3] - data2Re[3]);
455 data2Im[3] = -OnetoSqrt2 * (data2Re[3] + data2Im[3]);
456 data2Re[3] = tem;
457
458 temRe[0] = data1Re[0] + data2Re[0];
459 temRe[4] = data1Re[0] - data2Re[0];
460 temRe[1] = data1Re[1] + data2Re[1];
461 temRe[5] = data1Re[1] - data2Re[1];
462 temRe[2] = data1Re[2] + data2Re[2];
463 temRe[6] = data1Re[2] - data2Re[2];
464 temRe[3] = data1Re[3] + data2Re[3];
465 temRe[7] = data1Re[3] - data2Re[3];
466
467 temIm[0] = data1Im[0] + data2Im[0];
468 temIm[4] = data1Im[0] - data2Im[0];
469 temIm[1] = data1Im[1] + data2Im[1];
470 temIm[5] = data1Im[1] - data2Im[1];
471 temIm[2] = data1Im[2] + data2Im[2];
472 temIm[6] = data1Im[2] - data2Im[2];
473 temIm[3] = data1Im[3] + data2Im[3];
474 temIm[7] = data1Im[3] - data2Im[3];
475 } // End of function fft8().
476
477 private void fft10() {
478 float data1Re[] = new float[5];
479 float data1Im[] = new float[5];
480 float data2Re[] = new float[5];
481 float data2Im[] = new float[5];
482
483 // To improve the speed, use direct assaignment instead for loop here.
484 data1Re[0] = temRe[0];
485 data2Re[0] = temRe[5];
486 data1Re[1] = temRe[2];
487 data2Re[1] = temRe[7];
488 data1Re[2] = temRe[4];
489 data2Re[2] = temRe[9];
490 data1Re[3] = temRe[6];
491 data2Re[3] = temRe[1];
492 data1Re[4] = temRe[8];
493 data2Re[4] = temRe[3];
494
495 data1Im[0] = temIm[0];
496 data2Im[0] = temIm[5];
497 data1Im[1] = temIm[2];
498 data2Im[1] = temIm[7];
499 data1Im[2] = temIm[4];
500 data2Im[2] = temIm[9];
501 data1Im[3] = temIm[6];
502 data2Im[3] = temIm[1];
503 data1Im[4] = temIm[8];
504 data2Im[4] = temIm[3];
505
506 fft5(data1Re, data1Im);
507 fft5(data2Re, data2Im);
508
509 temRe[0] = data1Re[0] + data2Re[0];
510 temRe[5] = data1Re[0] - data2Re[0];
511 temRe[6] = data1Re[1] + data2Re[1];
512 temRe[1] = data1Re[1] - data2Re[1];
513 temRe[2] = data1Re[2] + data2Re[2];
514 temRe[7] = data1Re[2] - data2Re[2];
515 temRe[8] = data1Re[3] + data2Re[3];
516 temRe[3] = data1Re[3] - data2Re[3];
517 temRe[4] = data1Re[4] + data2Re[4];
518 temRe[9] = data1Re[4] - data2Re[4];
519
520 temIm[0] = data1Im[0] + data2Im[0];
521 temIm[5] = data1Im[0] - data2Im[0];
522 temIm[6] = data1Im[1] + data2Im[1];
523 temIm[1] = data1Im[1] - data2Im[1];
524 temIm[2] = data1Im[2] + data2Im[2];
525 temIm[7] = data1Im[2] - data2Im[2];
526 temIm[8] = data1Im[3] + data2Im[3];
527 temIm[3] = data1Im[3] - data2Im[3];
528 temIm[4] = data1Im[4] + data2Im[4];
529 temIm[9] = data1Im[4] - data2Im[4];
530 } // End of function fft10().
531
532 public double sqrt(double d) {
533 return Math.sqrt(d);
534 }
535
536 private void fftPrime(int radix) {
537 // Initial WRe, WIm.
538 float W = 2 * (float) Math.PI / radix;
539 float cosW = (float) Math.cos(W);
540 float sinW = -(float) Math.sin(W);
541 float WRe[] = new float[radix];
542 float WIm[] = new float[radix];
543
544 WRe[0] = 1;
545 WIm[0] = 0;
546 WRe[1] = cosW;
547 WIm[1] = sinW;
548
549 for (int i = 2; i < radix; i++) {
550 WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
551 WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
552 }
553
554 // FFT of prime length data, using DFT, can be improved in the future.
555 float rere, reim, imre, imim;
556 int j, k;
557 int max = (radix + 1) / 2;
558
559 float tem1Re[] = new float[max];
560 float tem1Im[] = new float[max];
561 float tem2Re[] = new float[max];
562 float tem2Im[] = new float[max];
563
564 for (j = 1; j < max; j++) {
565 tem1Re[j] = temRe[j] + temRe[radix - j];
566 tem1Im[j] = temIm[j] - temIm[radix - j];
567 tem2Re[j] = temRe[j] - temRe[radix - j];
568 tem2Im[j] = temIm[j] + temIm[radix - j];
569 }
570
571 for (j = 1; j < max; j++) {
572 temRe[j] = temRe[0];
573 temIm[j] = temIm[0];
574 temRe[radix - j] = temRe[0];
575 temIm[radix - j] = temIm[0];
576 k = j;
577 for (int i = 1; i < max; i++) {
578 rere = WRe[k] * tem1Re[i];
579 imim = WIm[k] * tem1Im[i];
580 reim = WRe[k] * tem2Im[i];
581 imre = WIm[k] * tem2Re[i];
582
583 temRe[radix - j] += rere + imim;
584 temIm[radix - j] += reim - imre;
585 temRe[j] += rere - imim;
586 temIm[j] += reim + imre;
587
588 k = k + j;
589 if (k >= radix)
590 k = k - radix;
591 }
592 }
593 for (j = 1; j < max; j++) {
594 temRe[0] = temRe[0] + tem1Re[j];
595 temIm[0] = temIm[0] + tem2Im[j];
596 }
597 } // End of function fftPrime().
598
599 }
600