1 public class fft2d extends Thread {
2 //Title: 2-d mixed radix FFT.
4 //Copyright: Copyright (c) 1998
6 //Company: University of Wisconsin-Milwaukee.
8 // . Use fft1d to perform fft2d.
10 // Code borrowed from :Java Digital Signal Processing book by Lyon and Rao
12 public Matrix data1, data2;
15 // Constructor: 2-d FFT of Complex data.
16 public fft2d(Matrix data1, Matrix data2, int x0, int x1) {
26 barr = new Barrier("128.195.175.84");
27 double tempdataRe[][];
28 double tempdataIm[][];
29 int rowlength, columnlength;
32 // Calculate FFT for each row of the data.
35 columnlength = data1.N;
36 tempdataRe = data1.dataRe;
37 tempdataIm = data1.dataIm;
40 fft1 = new fft1d(columnlength);
41 fft2 = new fft1d(rowlength);
42 for (int i = x0; i < x1; i++) {
44 double inputRe[] = tempdataRe[i]; //local array
45 double inputIm[] = tempdataIm[i];
46 fft(fft1, inputRe, inputIm);
51 Barrier.enterBarrier(barr);
56 transpose(tempdataRe,tempdataIm, data2.dataRe,data2.dataIm, rowlength, columnlength);
61 Barrier.enterBarrier(barr);
63 // Calculate FFT for each column of the data.
64 double transtempRe[][];
65 double transtempIm[][];
67 transtempRe = data2.dataRe;
68 transtempIm = data2.dataIm;
69 for (int j = start; j < end; j++) {
71 double inputRe[] = transtempRe[j]; //local array
72 double inputIm[] = transtempIm[j];
73 fft(fft2, inputRe, inputIm);
78 public void transpose(double[][] tempdataRe, double[][] tempdataIm, double[][] outputRe,
79 double[][] outputIm, int rowlength, int columnlength) {
80 for(int i = 0; i<rowlength; i++) {
81 double tRe[] = tempdataRe[i];
82 double tIm[] = tempdataIm[i];
83 for(int j = 0; j<columnlength; j++) {
84 outputRe[j][i] = tRe[j];
85 outputIm[j][i] = tIm[j];
90 public static void main(String[] args) {
95 NUM_THREADS=Integer.parseInt(args[0]);
97 SIZE = Integer.parseInt(args[1]);
100 System.printString("Num threads = " + NUM_THREADS + " SIZE= " + SIZE + "\n");
102 int[] mid = new int[8];
103 mid[0] = (128<<24)|(195<<16)|(136<<8)|162; //dw-10
104 mid[1] = (128<<24)|(195<<16)|(136<<8)|163; //dw-11
105 mid[2] = (128<<24)|(195<<16)|(136<<8)|164; //dw-12
106 mid[3] = (128<<24)|(195<<16)|(136<<8)|165; //dw-13
107 mid[4] = (128<<24)|(195<<16)|(136<<8)|166; //dw-14
108 mid[5] = (128<<24)|(195<<16)|(136<<8)|167; //dw-15
109 mid[6] = (128<<24)|(195<<16)|(136<<8)|168; //dw-16
110 mid[7] = (128<<24)|(195<<16)|(136<<8)|169; //dw-17
113 // Start Barrier Server
114 BarrierServer mybarr;
116 mybarr = global new BarrierServer(NUM_THREADS);
118 mybarr.start(mid[0]);
123 // Create threads to do FFT
126 // Set up data for FFT transform
127 data1 = global new Matrix(SIZE, SIZE);
128 data2 = global new Matrix(SIZE, SIZE);
129 data1.setValues(); //Input Matrix
130 data2.setZeros(); //Transpose Matrix
131 myfft2d = global new fft2d[NUM_THREADS];
132 int increment = SIZE/NUM_THREADS;
134 for(int i =0 ; i<NUM_THREADS; i++) {
135 if((i+1)==NUM_THREADS)
136 myfft2d[i] = global new fft2d(data1, data2, base, SIZE);
138 myfft2d[i] = global new fft2d(data1, data2, base, base+increment);
143 boolean waitfordone=true;
152 //Start a thread to compute each c[l,n]
153 for(int i = 0; i<NUM_THREADS; i++) {
160 //Wait for thread to finish
161 for(int i = 0; i<NUM_THREADS; i++) {
168 System.printString("2DFFT done! \n");
171 public static void fft(fft1d myfft, double inputRe[], double inputIm[]) {
173 double outputRe[] = myfft.outputRe;
174 double outputIm[] = myfft.outputIm;
175 // intermediate results
176 double temRe[] = myfft.temRe;
177 double temIm[] = myfft.temIm;
178 //Permute() operation
179 permute(myfft, outputRe, outputIm, inputRe, inputIm);
181 //System.printString("ready to twiddle");
182 for (int factorIndex = 0; factorIndex < myfft.NumofFactors; factorIndex++)
183 twiddle(factorIndex, myfft, temRe, temIm, outputRe, outputIm);
185 //System.printString("ready to copy");
186 // Copy the output[] data to input[], so the output can be
187 // returned in the input array.
188 for (int i = 0; i < myfft.N; i++) {
189 inputRe[i] = outputRe[i];
190 inputIm[i] = outputIm[i];
194 private static void permute(fft1d myfft, double[] outputRe, double[] outputIm, double[] inputRe, double[] inputIm) {
195 int count[] = new int[myfft.MaxFactorsNumber];
199 for (int i = 0; i < myfft.N - 1; i++) {
200 outputRe[i] = inputRe[k];
201 outputIm[i] = inputIm[k];
203 k = k + myfft.remain[j];
204 count[0] = count[0] + 1;
205 while (count[j] >= myfft.factors[j]) {
211 tmp = myfft.remain[j - 1];
212 k = k - tmp + myfft.remain[j + 1];
214 count[j] = count[j] + 1;
217 outputRe[myfft.N - 1] = inputRe[myfft.N - 1];
218 outputIm[myfft.N - 1] = inputIm[myfft.N - 1];
219 } // End of function permute().
221 private static void twiddle(int factorIndex, fft1d myfft, double[] temRe, double[] temIm,
222 double[] outputRe, double[] outputIm) {
224 int sofarRadix = myfft.sofar[factorIndex];
225 int radix = myfft.factors[factorIndex];
226 int remainRadix = myfft.remain[factorIndex];
228 double tem; // Temporary variable to do data exchange.
230 double W = 2 * (double) Math.setPI() / (sofarRadix * radix);
231 double cosW = (double) Math.cos(W);
232 double sinW = -(double) Math.sin(W);
234 double twiddleRe[] = new double[radix];
235 double twiddleIm[] = new double[radix];
236 double twRe = 1.0f, twIm = 0f;
238 //Initialize twiddle addBk.address variables.
239 int dataOffset = 0, groupOffset = 0, address = 0;
241 for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
242 //System.printString("datano="+dataNo);
243 if (sofarRadix > 1) {
248 for (int i = 2; i < radix; i++) {
249 twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
250 twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
252 tem = cosW * twRe - sinW * twIm;
253 twIm = sinW * twRe + cosW * twIm;
256 for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
257 //System.printString("groupNo="+groupNo);
258 if ((sofarRadix > 1) && (dataNo > 0)) {
259 temRe[0] = outputRe[address];
260 temIm[0] = outputIm[address];
263 address = address + sofarRadix;
264 temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
265 twiddleIm[blockIndex] * outputIm[address];
266 temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
267 twiddleIm[blockIndex] * outputRe[address];
269 } while (blockIndex < radix);
271 for (int i = 0; i < radix; i++) {
272 //System.printString("temRe.length="+temRe.length);
273 //System.printString("i = "+i);
274 temRe[i] = outputRe[address];
275 temIm[i] = outputIm[address];
276 address += sofarRadix;
279 //System.printString("radix="+radix);
281 tem = temRe[0] + temRe[1];
282 temRe[1] = temRe[0] - temRe[1];
284 tem = temIm[0] + temIm[1];
285 temIm[1] = temIm[0] - temIm[1];
287 } else if( radix == 3) {
288 double t1Re = temRe[1] + temRe[2];
289 double t1Im = temIm[1] + temIm[2];
290 temRe[0] = temRe[0] + t1Re;
291 temIm[0] = temIm[0] + t1Im;
293 double m1Re = myfft.cos2to3PI * t1Re;
294 double m1Im = myfft.cos2to3PI * t1Im;
295 double m2Re = myfft.sin2to3PI * (temIm[1] - temIm[2]);
296 double m2Im = myfft.sin2to3PI * (temRe[2] - temRe[1]);
297 double s1Re = temRe[0] + m1Re;
298 double s1Im = temIm[0] + m1Im;
300 temRe[1] = s1Re + m2Re;
301 temIm[1] = s1Im + m2Im;
302 temRe[2] = s1Re - m2Re;
303 temIm[2] = s1Im - m2Im;
304 } else if(radix == 4) {
306 } else if(radix == 5) {
307 fft5(myfft, temRe, temIm);
308 } else if(radix == 8) {
309 fft8(myfft, temRe, temIm);
310 } else if(radix == 10) {
311 fft10(myfft, temRe, temIm);
313 fftPrime(radix, temRe, temIm);
315 address = groupOffset;
316 for (int i = 0; i < radix; i++) {
317 outputRe[address] = temRe[i];
318 outputIm[address] = temIm[i];
319 address += sofarRadix;
321 groupOffset += sofarRadix * radix;
322 address = groupOffset;
324 groupOffset = ++dataOffset;
325 address = groupOffset;
327 } //twiddle operation
329 // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
330 private static void fft4(double dataRe[], double dataIm[]) {
331 double t1Re,t1Im, t2Re,t2Im;
332 double m2Re,m2Im, m3Re,m3Im;
334 t1Re = dataRe[0] + dataRe[2];
335 t1Im = dataIm[0] + dataIm[2];
336 t2Re = dataRe[1] + dataRe[3];
337 t2Im = dataIm[1] + dataIm[3];
339 m2Re = dataRe[0] - dataRe[2];
340 m2Im = dataIm[0] - dataIm[2];
341 m3Re = dataIm[1] - dataIm[3];
342 m3Im = dataRe[3] - dataRe[1];
344 dataRe[0] = t1Re + t2Re;
345 dataIm[0] = t1Im + t2Im;
346 dataRe[2] = t1Re - t2Re;
347 dataIm[2] = t1Im - t2Im;
348 dataRe[1] = m2Re + m3Re;
349 dataIm[1] = m2Im + m3Im;
350 dataRe[3] = m2Re - m3Re;
351 dataIm[3] = m2Im - m3Im;
352 } // End of function fft4().
354 // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
355 private static void fft5(fft1d myfft, double dataRe[], double dataIm[]) {
356 double t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
357 double m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
358 double s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
360 t1Re = dataRe[1] + dataRe[4];
361 t1Im = dataIm[1] + dataIm[4];
362 t2Re = dataRe[2] + dataRe[3];
363 t2Im = dataIm[2] + dataIm[3];
364 t3Re = dataRe[1] - dataRe[4];
365 t3Im = dataIm[1] - dataIm[4];
366 t4Re = dataRe[3] - dataRe[2];
367 t4Im = dataIm[3] - dataIm[2];
371 dataRe[0] = dataRe[0] + t5Re;
372 dataIm[0] = dataIm[0] + t5Im;
374 m1Re = myfft.c51 * t5Re;
375 m1Im = myfft.c51 * t5Im;
376 m2Re = myfft.c52 * (t1Re - t2Re);
377 m2Im = myfft.c52 * (t1Im - t2Im);
378 m3Re = -(myfft.c53) * (t3Im + t4Im);
379 m3Im = myfft.c53 * (t3Re + t4Re);
380 m4Re = -(myfft.c54) * t4Im;
381 m4Im = myfft.c54 * t4Re;
382 m5Re = -(myfft.c55) * t3Im;
383 m5Im = myfft.c55 * t3Re;
389 s1Re = dataRe[0] + m1Re;
390 s1Im = dataIm[0] + m1Im;
396 dataRe[1] = s2Re + s3Re;
397 dataIm[1] = s2Im + s3Im;
398 dataRe[2] = s4Re + s5Re;
399 dataIm[2] = s4Im + s5Im;
400 dataRe[3] = s4Re - s5Re;
401 dataIm[3] = s4Im - s5Im;
402 dataRe[4] = s2Re - s3Re;
403 dataIm[4] = s2Im - s3Im;
404 } // End of function fft5().
406 private static void fft8(fft1d myfft, double[] temRe, double[] temIm) {
407 double data1Re[] = new double[4];
408 double data1Im[] = new double[4];
409 double data2Re[] = new double[4];
410 double data2Im[] = new double[4];
413 // To improve the speed, use direct assaignment instead for loop here.
414 data1Re[0] = temRe[0];
415 data2Re[0] = temRe[1];
416 data1Re[1] = temRe[2];
417 data2Re[1] = temRe[3];
418 data1Re[2] = temRe[4];
419 data2Re[2] = temRe[5];
420 data1Re[3] = temRe[6];
421 data2Re[3] = temRe[7];
423 data1Im[0] = temIm[0];
424 data2Im[0] = temIm[1];
425 data1Im[1] = temIm[2];
426 data2Im[1] = temIm[3];
427 data1Im[2] = temIm[4];
428 data2Im[2] = temIm[5];
429 data1Im[3] = temIm[6];
430 data2Im[3] = temIm[7];
432 fft4(data1Re, data1Im);
433 fft4(data2Re, data2Im);
435 tem = myfft.OnetoSqrt2 * (data2Re[1] + data2Im[1]);
436 data2Im[1] = myfft.OnetoSqrt2 * (data2Im[1] - data2Re[1]);
439 data2Im[2] = -data2Re[2];
441 tem = myfft.OnetoSqrt2 * (data2Im[3] - data2Re[3]);
442 data2Im[3] = -(myfft.OnetoSqrt2) * (data2Re[3] + data2Im[3]);
445 temRe[0] = data1Re[0] + data2Re[0];
446 temRe[4] = data1Re[0] - data2Re[0];
447 temRe[1] = data1Re[1] + data2Re[1];
448 temRe[5] = data1Re[1] - data2Re[1];
449 temRe[2] = data1Re[2] + data2Re[2];
450 temRe[6] = data1Re[2] - data2Re[2];
451 temRe[3] = data1Re[3] + data2Re[3];
452 temRe[7] = data1Re[3] - data2Re[3];
454 temIm[0] = data1Im[0] + data2Im[0];
455 temIm[4] = data1Im[0] - data2Im[0];
456 temIm[1] = data1Im[1] + data2Im[1];
457 temIm[5] = data1Im[1] - data2Im[1];
458 temIm[2] = data1Im[2] + data2Im[2];
459 temIm[6] = data1Im[2] - data2Im[2];
460 temIm[3] = data1Im[3] + data2Im[3];
461 temIm[7] = data1Im[3] - data2Im[3];
462 } // End of function fft8().
464 private static void fft10(fft1d myfft, double[] temRe, double[] temIm) {
465 double data1Re[] = new double[5];
466 double data1Im[] = new double[5];
467 double data2Re[] = new double[5];
468 double data2Im[] = new double[5];
470 // To improve the speed, use direct assaignment instead for loop here.
471 data1Re[0] = temRe[0];
472 data2Re[0] = temRe[5];
473 data1Re[1] = temRe[2];
474 data2Re[1] = temRe[7];
475 data1Re[2] = temRe[4];
476 data2Re[2] = temRe[9];
477 data1Re[3] = temRe[6];
478 data2Re[3] = temRe[1];
479 data1Re[4] = temRe[8];
480 data2Re[4] = temRe[3];
482 data1Im[0] = temIm[0];
483 data2Im[0] = temIm[5];
484 data1Im[1] = temIm[2];
485 data2Im[1] = temIm[7];
486 data1Im[2] = temIm[4];
487 data2Im[2] = temIm[9];
488 data1Im[3] = temIm[6];
489 data2Im[3] = temIm[1];
490 data1Im[4] = temIm[8];
491 data2Im[4] = temIm[3];
493 fft5(myfft, data1Re, data1Im);
494 fft5(myfft, data2Re, data2Im);
496 temRe[0] = data1Re[0] + data2Re[0];
497 temRe[5] = data1Re[0] - data2Re[0];
498 temRe[6] = data1Re[1] + data2Re[1];
499 temRe[1] = data1Re[1] - data2Re[1];
500 temRe[2] = data1Re[2] + data2Re[2];
501 temRe[7] = data1Re[2] - data2Re[2];
502 temRe[8] = data1Re[3] + data2Re[3];
503 temRe[3] = data1Re[3] - data2Re[3];
504 temRe[4] = data1Re[4] + data2Re[4];
505 temRe[9] = data1Re[4] - data2Re[4];
507 temIm[0] = data1Im[0] + data2Im[0];
508 temIm[5] = data1Im[0] - data2Im[0];
509 temIm[6] = data1Im[1] + data2Im[1];
510 temIm[1] = data1Im[1] - data2Im[1];
511 temIm[2] = data1Im[2] + data2Im[2];
512 temIm[7] = data1Im[2] - data2Im[2];
513 temIm[8] = data1Im[3] + data2Im[3];
514 temIm[3] = data1Im[3] - data2Im[3];
515 temIm[4] = data1Im[4] + data2Im[4];
516 temIm[9] = data1Im[4] - data2Im[4];
517 } // End of function fft10().
519 private static void fftPrime(int radix, double[] temRe, double[] temIm) {
521 double W = 2 * (double) Math.setPI() / radix;
522 double cosW = (double) Math.cos(W);
523 double sinW = -(double) Math.sin(W);
524 double WRe[] = new double[radix];
525 double WIm[] = new double[radix];
532 for (int i = 2; i < radix; i++) {
533 WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
534 WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
537 // FFT of prime length data, using DFT, can be improved in the future.
538 double rere, reim, imre, imim;
540 int max = (radix + 1) / 2;
542 double tem1Re[] = new double[max];
543 double tem1Im[] = new double[max];
544 double tem2Re[] = new double[max];
545 double tem2Im[] = new double[max];
547 for (j = 1; j < max; j++) {
548 tem1Re[j] = temRe[j] + temRe[radix - j];
549 tem1Im[j] = temIm[j] - temIm[radix - j];
550 tem2Re[j] = temRe[j] - temRe[radix - j];
551 tem2Im[j] = temIm[j] + temIm[radix - j];
554 for (j = 1; j < max; j++) {
557 temRe[radix - j] = temRe[0];
558 temIm[radix - j] = temIm[0];
560 for (int i = 1; i < max; i++) {
561 rere = WRe[k] * tem1Re[i];
562 imim = WIm[k] * tem1Im[i];
563 reim = WRe[k] * tem2Im[i];
564 imre = WIm[k] * tem2Re[i];
566 temRe[radix - j] += rere + imim;
567 temIm[radix - j] += reim - imre;
568 temRe[j] += rere - imim;
569 temIm[j] += reim + imre;
576 for (j = 1; j < max; j++) {
577 temRe[0] = temRe[0] + tem1Re[j];
578 temIm[0] = temIm[0] + tem2Im[j];
580 } // End of function fftPrime().