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 fft1d fft1, fft2;
14 public int x0, x1, y0, y1;
15 public double inputRe[], inputIm[];
17 // Constructor: 2-d FFT of Complex data.
18 public fft2d(double[] inputRe, double[] inputIm, Matrix data, fft1d fft1, fft1d fft2, int x0, int x1, int y0, int y1) {
26 this.inputRe = inputRe;
27 this.inputIm = inputIm;
32 barr = new Barrier("128.195.175.84");
33 double tempdataRe[][];
34 double tempdataIm[][];
37 int rowlength, columnlength;
40 rowlength = data.M; //height
41 columnlength = data.N; //width
42 tempdataRe = data.dataRe;
43 tempdataIm = data.dataIm;
45 // Calculate FFT for each row of the data.
46 //System.printString("x0= " + x0 + " x1= " + x1 + " y0= "+ y0 + " y1= " + y1 + " width = " + columnlength + " height= " + rowlength+ "\n");
47 for (int i = x0; i < x1; i++) {
49 if(columnlength != N) {
50 System.printString("Error: the length of real part & imaginary part " + "of the input to 1-d FFT are different");
53 //Permute() operation on fft1
55 double inputRe[] = tempdataRe[i]; //local array
56 double inputIm[] = tempdataIm[i];
58 double outputRe[] = fft1.outputRe; //local array
59 double outputIm[] = fft1.outputIm;
60 double temRe[] = fft1.temRe; // intermediate results
61 double temIm[] = fft1.temIm;
62 int count[] = new int[fft1.MaxFactorsNumber];
65 for(int a = 0; a < N-1; a++) {
66 outputRe[a] = inputRe[k];
67 outputIm[a] = inputIm[k];
69 k = k + fft1.remain[j];
70 count[0] = count[0] + 1;
71 while (count[j] >= fft1.factors[j]) {
77 tmp = fft1.remain[j - 1];
78 k = k - tmp + fft1.remain[j + 1];
80 count[j] = count[j] + 1;
83 outputRe[N - 1] = inputRe[N - 1];
84 outputIm[N - 1] = inputIm[N - 1];
86 //Twiddle oeration on fft1
87 for (int factorIndex = 0; factorIndex < fft1.NumofFactors; factorIndex++) {
88 twiddle(factorIndex, fft1, temRe, temIm, outputRe, outputIm);
90 // Copy the output[] data to input[], so the output can be
91 // returned in the input array.
92 for (int b = 0; b < N; b++) {
93 inputRe[b] = outputRe[b];
94 inputIm[b] = outputIm[b];
101 Barrier.enterBarrier(barr);
105 mytemRe = new double[columnlength][rowlength];
106 mytemIm = new double[columnlength][rowlength];
107 for(int i = x0; i<x1; i++) {
108 double tRe[] = tempdataRe[i];
109 double tIm[] = tempdataIm[i];
110 for(int j = y0; j<y1; j++) {
111 mytemRe[j][i] = tRe[j];
112 mytemIm[j][i] = tIm[j];
118 Barrier.enterBarrier(barr);
120 // Calculate FFT for each column of the data.
122 for (int j = y0; j < y1; j++) {
125 System.printString("Error: the length of real part & imaginary part " + "of the input to 1-d FFT are different");
128 //Permute() operation on fft2
130 double inputRe[] = mytemRe[j]; //local array
131 double inputIm[] = mytemIm[j];
133 double outputRe[] = fft2.outputRe; //local array
134 double outputIm[] = fft2.outputIm;
135 double temRe[] = fft2.temRe; // intermediate results
136 double temIm[] = fft2.temIm;
137 int count[] = new int[fft2.MaxFactorsNumber];
140 for(int a = 0; a < N-1; a++) {
141 outputRe[a] = inputRe[k];
142 outputIm[a] = inputIm[k];
144 k = k + fft2.remain[r];
145 count[0] = count[0] + 1;
146 while (count[r] >= fft2.factors[r]) {
152 tmp = fft2.remain[r - 1];
153 k = k - tmp + fft2.remain[r + 1];
155 count[r] = count[r] + 1;
158 outputRe[N - 1] = inputRe[N - 1];
159 outputIm[N - 1] = inputIm[N - 1];
161 //Twiddle oeration on fft2
162 for (int factorIndex = 0; factorIndex < fft2.NumofFactors; factorIndex++) {
163 twiddle(factorIndex, fft2, temRe, temIm, outputRe, outputIm);
165 // Copy the output[] data to input[], so the output can be
166 // returned in the input array.
167 for (int b = 0; b < N; b++) {
168 inputRe[b] = outputRe[b];
169 inputIm[b] = outputIm[b];
176 Barrier.enterBarrier(barr);
179 // Copy the result to input[], so the output can be
180 // returned in the input array.
182 for (int j = y0; j < y1; j++) {
183 double tRe[] = mytemRe[j];
184 double tIm[] = mytemIm[j];
185 for (int i = x0; i < x1; i++) {
186 inputRe[i* data.N + j] = tRe[i];
187 inputIm[i* data.N + j] = tIm[i];
195 //("ready to twiddle");
196 private void twiddle(int factorIndex, fft1d myfft, double[] temRe, double[] temIm,
197 double[] outputRe, double[] outputIm) {
199 int sofarRadix = myfft.sofar[factorIndex];
200 int radix = myfft.factors[factorIndex];
201 int remainRadix = myfft.remain[factorIndex];
203 double tem; // Temporary variable to do data exchange.
205 double W = 2 * (double) Math.setPI() / (sofarRadix * radix);
206 double cosW = (double) Math.cos(W);
207 double sinW = -(double) Math.sin(W);
209 double twiddleRe[] = new double[radix];
210 double twiddleIm[] = new double[radix];
211 double twRe = 1.0f, twIm = 0f;
213 //Initialize twiddle addBk.address variables.
214 int dataOffset = 0, groupOffset = 0, address = 0;
216 for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
217 //System.printString("datano="+dataNo);
218 if (sofarRadix > 1) {
223 for (int i = 2; i < radix; i++) {
224 twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
225 twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
227 tem = cosW * twRe - sinW * twIm;
228 twIm = sinW * twRe + cosW * twIm;
231 for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
232 //System.printString("groupNo="+groupNo);
233 if ((sofarRadix > 1) && (dataNo > 0)) {
234 temRe[0] = outputRe[address];
235 temIm[0] = outputIm[address];
238 address = address + sofarRadix;
239 temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
240 twiddleIm[blockIndex] * outputIm[address];
241 temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
242 twiddleIm[blockIndex] * outputRe[address];
244 } while (blockIndex < radix);
246 for (int i = 0; i < radix; i++) {
247 //System.printString("temRe.length="+temRe.length);
248 //System.printString("i = "+i);
249 temRe[i] = outputRe[address];
250 temIm[i] = outputIm[address];
251 address += sofarRadix;
253 //System.printString("radix="+radix);
255 tem = temRe[0] + temRe[1];
256 temRe[1] = temRe[0] - temRe[1];
258 tem = temIm[0] + temIm[1];
259 temIm[1] = temIm[0] - temIm[1];
261 } else if( radix == 3) {
262 double t1Re = temRe[1] + temRe[2];
263 double t1Im = temIm[1] + temIm[2];
264 temRe[0] = temRe[0] + t1Re;
265 temIm[0] = temIm[0] + t1Im;
267 double m1Re = myfft.cos2to3PI * t1Re;
268 double m1Im = myfft.cos2to3PI * t1Im;
269 double m2Re = myfft.sin2to3PI * (temIm[1] - temIm[2]);
270 double m2Im = myfft.sin2to3PI * (temRe[2] - temRe[1]);
271 double s1Re = temRe[0] + m1Re;
272 double s1Im = temIm[0] + m1Im;
274 temRe[1] = s1Re + m2Re;
275 temIm[1] = s1Im + m2Im;
276 temRe[2] = s1Re - m2Re;
277 temIm[2] = s1Im - m2Im;
278 } else if(radix == 4) {
280 } else if(radix == 5) {
281 fft5(myfft, temRe, temIm);
282 } else if(radix == 8) {
283 fft8(myfft, temRe, temIm);
284 } else if(radix == 10) {
285 fft10(myfft, temRe, temIm);
287 fftPrime(radix, temRe, temIm);
289 address = groupOffset;
290 for (int i = 0; i < radix; i++) {
291 outputRe[address] = temRe[i];
292 outputIm[address] = temIm[i];
293 address += sofarRadix;
295 groupOffset += sofarRadix * radix;
296 address = groupOffset;
298 groupOffset = ++dataOffset;
299 address = groupOffset;
301 } //twiddle operation
303 // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
304 private void fft4(double dataRe[], double dataIm[]) {
305 double t1Re,t1Im, t2Re,t2Im;
306 double m2Re,m2Im, m3Re,m3Im;
308 t1Re = dataRe[0] + dataRe[2];
309 t1Im = dataIm[0] + dataIm[2];
310 t2Re = dataRe[1] + dataRe[3];
311 t2Im = dataIm[1] + dataIm[3];
313 m2Re = dataRe[0] - dataRe[2];
314 m2Im = dataIm[0] - dataIm[2];
315 m3Re = dataIm[1] - dataIm[3];
316 m3Im = dataRe[3] - dataRe[1];
318 dataRe[0] = t1Re + t2Re;
319 dataIm[0] = t1Im + t2Im;
320 dataRe[2] = t1Re - t2Re;
321 dataIm[2] = t1Im - t2Im;
322 dataRe[1] = m2Re + m3Re;
323 dataIm[1] = m2Im + m3Im;
324 dataRe[3] = m2Re - m3Re;
325 dataIm[3] = m2Im - m3Im;
326 } // End of function fft4().
328 // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
329 private void fft5(fft1d myfft, double dataRe[], double dataIm[]) {
330 double t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
331 double m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
332 double s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
334 t1Re = dataRe[1] + dataRe[4];
335 t1Im = dataIm[1] + dataIm[4];
336 t2Re = dataRe[2] + dataRe[3];
337 t2Im = dataIm[2] + dataIm[3];
338 t3Re = dataRe[1] - dataRe[4];
339 t3Im = dataIm[1] - dataIm[4];
340 t4Re = dataRe[3] - dataRe[2];
341 t4Im = dataIm[3] - dataIm[2];
345 dataRe[0] = dataRe[0] + t5Re;
346 dataIm[0] = dataIm[0] + t5Im;
348 m1Re = myfft.c51 * t5Re;
349 m1Im = myfft.c51 * t5Im;
350 m2Re = myfft.c52 * (t1Re - t2Re);
351 m2Im = myfft.c52 * (t1Im - t2Im);
352 m3Re = -(myfft.c53) * (t3Im + t4Im);
353 m3Im = myfft.c53 * (t3Re + t4Re);
354 m4Re = -(myfft.c54) * t4Im;
355 m4Im = myfft.c54 * t4Re;
356 m5Re = -(myfft.c55) * t3Im;
357 m5Im = myfft.c55 * t3Re;
363 s1Re = dataRe[0] + m1Re;
364 s1Im = dataIm[0] + m1Im;
370 dataRe[1] = s2Re + s3Re;
371 dataIm[1] = s2Im + s3Im;
372 dataRe[2] = s4Re + s5Re;
373 dataIm[2] = s4Im + s5Im;
374 dataRe[3] = s4Re - s5Re;
375 dataIm[3] = s4Im - s5Im;
376 dataRe[4] = s2Re - s3Re;
377 dataIm[4] = s2Im - s3Im;
378 } // End of function fft5().
380 private void fft8(fft1d myfft, double[] temRe, double[] temIm) {
381 double data1Re[] = new double[4];
382 double data1Im[] = new double[4];
383 double data2Re[] = new double[4];
384 double data2Im[] = new double[4];
387 // To improve the speed, use direct assaignment instead for loop here.
388 data1Re[0] = temRe[0];
389 data2Re[0] = temRe[1];
390 data1Re[1] = temRe[2];
391 data2Re[1] = temRe[3];
392 data1Re[2] = temRe[4];
393 data2Re[2] = temRe[5];
394 data1Re[3] = temRe[6];
395 data2Re[3] = temRe[7];
397 data1Im[0] = temIm[0];
398 data2Im[0] = temIm[1];
399 data1Im[1] = temIm[2];
400 data2Im[1] = temIm[3];
401 data1Im[2] = temIm[4];
402 data2Im[2] = temIm[5];
403 data1Im[3] = temIm[6];
404 data2Im[3] = temIm[7];
406 fft4(data1Re, data1Im);
407 fft4(data2Re, data2Im);
409 tem = myfft.OnetoSqrt2 * (data2Re[1] + data2Im[1]);
410 data2Im[1] = myfft.OnetoSqrt2 * (data2Im[1] - data2Re[1]);
413 data2Im[2] = -data2Re[2];
415 tem = myfft.OnetoSqrt2 * (data2Im[3] - data2Re[3]);
416 data2Im[3] = -(myfft.OnetoSqrt2) * (data2Re[3] + data2Im[3]);
419 temRe[0] = data1Re[0] + data2Re[0];
420 temRe[4] = data1Re[0] - data2Re[0];
421 temRe[1] = data1Re[1] + data2Re[1];
422 temRe[5] = data1Re[1] - data2Re[1];
423 temRe[2] = data1Re[2] + data2Re[2];
424 temRe[6] = data1Re[2] - data2Re[2];
425 temRe[3] = data1Re[3] + data2Re[3];
426 temRe[7] = data1Re[3] - data2Re[3];
428 temIm[0] = data1Im[0] + data2Im[0];
429 temIm[4] = data1Im[0] - data2Im[0];
430 temIm[1] = data1Im[1] + data2Im[1];
431 temIm[5] = data1Im[1] - data2Im[1];
432 temIm[2] = data1Im[2] + data2Im[2];
433 temIm[6] = data1Im[2] - data2Im[2];
434 temIm[3] = data1Im[3] + data2Im[3];
435 temIm[7] = data1Im[3] - data2Im[3];
436 } // End of function fft8().
438 private void fft10(fft1d myfft, double[] temRe, double[] temIm) {
439 double data1Re[] = new double[5];
440 double data1Im[] = new double[5];
441 double data2Re[] = new double[5];
442 double data2Im[] = new double[5];
444 // To improve the speed, use direct assaignment instead for loop here.
445 data1Re[0] = temRe[0];
446 data2Re[0] = temRe[5];
447 data1Re[1] = temRe[2];
448 data2Re[1] = temRe[7];
449 data1Re[2] = temRe[4];
450 data2Re[2] = temRe[9];
451 data1Re[3] = temRe[6];
452 data2Re[3] = temRe[1];
453 data1Re[4] = temRe[8];
454 data2Re[4] = temRe[3];
456 data1Im[0] = temIm[0];
457 data2Im[0] = temIm[5];
458 data1Im[1] = temIm[2];
459 data2Im[1] = temIm[7];
460 data1Im[2] = temIm[4];
461 data2Im[2] = temIm[9];
462 data1Im[3] = temIm[6];
463 data2Im[3] = temIm[1];
464 data1Im[4] = temIm[8];
465 data2Im[4] = temIm[3];
467 fft5(myfft, data1Re, data1Im);
468 fft5(myfft, data2Re, data2Im);
470 temRe[0] = data1Re[0] + data2Re[0];
471 temRe[5] = data1Re[0] - data2Re[0];
472 temRe[6] = data1Re[1] + data2Re[1];
473 temRe[1] = data1Re[1] - data2Re[1];
474 temRe[2] = data1Re[2] + data2Re[2];
475 temRe[7] = data1Re[2] - data2Re[2];
476 temRe[8] = data1Re[3] + data2Re[3];
477 temRe[3] = data1Re[3] - data2Re[3];
478 temRe[4] = data1Re[4] + data2Re[4];
479 temRe[9] = data1Re[4] - data2Re[4];
481 temIm[0] = data1Im[0] + data2Im[0];
482 temIm[5] = data1Im[0] - data2Im[0];
483 temIm[6] = data1Im[1] + data2Im[1];
484 temIm[1] = data1Im[1] - data2Im[1];
485 temIm[2] = data1Im[2] + data2Im[2];
486 temIm[7] = data1Im[2] - data2Im[2];
487 temIm[8] = data1Im[3] + data2Im[3];
488 temIm[3] = data1Im[3] - data2Im[3];
489 temIm[4] = data1Im[4] + data2Im[4];
490 temIm[9] = data1Im[4] - data2Im[4];
491 } // End of function fft10().
493 private void fftPrime(int radix, double[] temRe, double[] temIm) {
495 double W = 2 * (double) Math.setPI() / radix;
496 double cosW = (double) Math.cos(W);
497 double sinW = -(double) Math.sin(W);
498 double WRe[] = new double[radix];
499 double WIm[] = new double[radix];
506 for (int i = 2; i < radix; i++) {
507 WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
508 WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
511 // FFT of prime length data, using DFT, can be improved in the future.
512 double rere, reim, imre, imim;
514 int max = (radix + 1) / 2;
516 double tem1Re[] = new double[max];
517 double tem1Im[] = new double[max];
518 double tem2Re[] = new double[max];
519 double tem2Im[] = new double[max];
521 for (j = 1; j < max; j++) {
522 tem1Re[j] = temRe[j] + temRe[radix - j];
523 tem1Im[j] = temIm[j] - temIm[radix - j];
524 tem2Re[j] = temRe[j] - temRe[radix - j];
525 tem2Im[j] = temIm[j] + temIm[radix - j];
528 for (j = 1; j < max; j++) {
531 temRe[radix - j] = temRe[0];
532 temIm[radix - j] = temIm[0];
534 for (int i = 1; i < max; i++) {
535 rere = WRe[k] * tem1Re[i];
536 imim = WIm[k] * tem1Im[i];
537 reim = WRe[k] * tem2Im[i];
538 imre = WIm[k] * tem2Re[i];
540 temRe[radix - j] += rere + imim;
541 temIm[radix - j] += reim - imre;
542 temRe[j] += rere - imim;
543 temIm[j] += reim + imre;
550 for (j = 1; j < max; j++) {
551 temRe[0] = temRe[0] + tem1Re[j];
552 temIm[0] = temIm[0] + tem2Im[j];
554 } // End of function fftPrime().
556 public static void main(String[] args) {
561 NUM_THREADS=Integer.parseInt(args[0]);
563 SIZE = Integer.parseInt(args[1]);
566 System.printString("Num threads = " + NUM_THREADS + " SIZE= " + SIZE + "\n");
569 // Matrix inputRe, inputIm;
574 inputRe = global new double[SIZE];
575 inputIm = global new double[SIZE];
577 for(int i = 0; i<SIZE; i++){
585 System.printString("Element 231567 is " + (int)inputRe[231567]+ "\n");
586 System.printString("Element 10 is " + (int)inputIm[10]+ "\n");
590 int[] mid = new int[8];
591 mid[0] = (128<<24)|(195<<16)|(175<<8)|84; //dw-10
592 mid[1] = (128<<24)|(195<<16)|(175<<8)|85; //dw-11
593 mid[2] = (128<<24)|(195<<16)|(175<<8)|86; //dw-12
594 mid[3] = (128<<24)|(195<<16)|(175<<8)|87; //dw-13
595 mid[4] = (128<<24)|(195<<16)|(175<<8)|88; //dw-14
596 mid[5] = (128<<24)|(195<<16)|(175<<8)|89; //dw-15
597 mid[6] = (128<<24)|(195<<16)|(175<<8)|90; //dw-16
598 mid[7] = (128<<24)|(195<<16)|(175<<8)|91; //dw-17
600 // Start Barrier Server
601 BarrierServer mybarr;
603 mybarr = global new BarrierServer(NUM_THREADS);
605 mybarr.start(mid[0]);
607 // Width and height of 2-d matrix inputRe or inputIm.
610 int Relength, Imlength;
612 height = inputRe.length / width;
613 Relength = inputRe.length;
614 Imlength = inputIm.length;
617 //System.printString("Initialized width and height\n");
620 // First make sure inputRe & inputIm are of the same length in terms of columns
621 if (Relength != Imlength) {
622 System.printString("Error: the length of real part & imaginary part " +
623 "of the input to 2-d FFT are different");
627 fft1 = global new fft1d(width);
628 fft2 = global new fft1d(height);
629 // Set up data for FFT transform
630 data = global new Matrix(height, width);
631 data.setValues(inputRe, inputIm);
634 // Create threads to do FFT
637 myfft2d = global new fft2d[NUM_THREADS];
638 int increment = height/NUM_THREADS;
640 for(int i =0 ; i<NUM_THREADS; i++) {
641 if((i+1)==NUM_THREADS)
642 myfft2d[i] = global new fft2d(inputRe, inputIm, data, fft1, fft2, base, increment, 0, width);
644 myfft2d[i] = global new fft2d(inputRe, inputIm, data, fft1, fft2, base, base+increment, 0, width);
649 boolean waitfordone=true;
658 //Start a thread to compute each c[l,n]
659 for(int i = 0; i<NUM_THREADS; i++) {
666 //Wait for thread to finish
667 for(int i = 0; i<NUM_THREADS; i++) {
675 System.printString("2DFFT done! \n");
678 System.printString("Element 231567 is " + (int)inputRe[231567]+ "\n");
679 System.printString("Element 10 is " + (int)inputIm[10]+ "\n");
685 // Copy the result to input[], so the output can be
686 // returned in the input array.
689 for (int i = 0; i < height; i++) {
690 for (int j = 0; j < width; j++) {
691 System.printString((int)inputRe[i * width + j]+ "\n");
692 System.printString((int)inputIm[i * width + j]+ "\n");