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)|(175<<8)|84; //dw-10
104 mid[1] = (128<<24)|(195<<16)|(175<<8)|85; //dw-11
105 mid[2] = (128<<24)|(195<<16)|(175<<8)|86; //dw-12
106 mid[3] = (128<<24)|(195<<16)|(175<<8)|87; //dw-13
107 mid[4] = (128<<24)|(195<<16)|(175<<8)|88; //dw-14
108 mid[5] = (128<<24)|(195<<16)|(175<<8)|89; //dw-15
109 mid[6] = (128<<24)|(195<<16)|(175<<8)|90; //dw-16
110 mid[7] = (128<<24)|(195<<16)|(175<<8)|91; //dw-17
112 // Start Barrier Server
113 BarrierServer mybarr;
115 mybarr = global new BarrierServer(NUM_THREADS);
117 mybarr.start(mid[0]);
122 // Create threads to do FFT
125 // Set up data for FFT transform
126 data1 = global new Matrix(SIZE, SIZE);
127 data2 = global new Matrix(SIZE, SIZE);
128 data1.setValues(); //Input Matrix
129 data2.setZeros(); //Transpose Matrix
130 myfft2d = global new fft2d[NUM_THREADS];
131 int increment = SIZE/NUM_THREADS;
133 for(int i =0 ; i<NUM_THREADS; i++) {
134 if((i+1)==NUM_THREADS)
135 myfft2d[i] = global new fft2d(data1, data2, base, SIZE);
137 myfft2d[i] = global new fft2d(data1, data2, base, base+increment);
142 boolean waitfordone=true;
151 //Start a thread to compute each c[l,n]
152 for(int i = 0; i<NUM_THREADS; i++) {
159 //Wait for thread to finish
160 for(int i = 0; i<NUM_THREADS; i++) {
167 System.printString("2DFFT done! \n");
170 public static void fft(fft1d myfft, double inputRe[], double inputIm[]) {
172 double outputRe[] = myfft.outputRe;
173 double outputIm[] = myfft.outputIm;
174 // intermediate results
175 double temRe[] = myfft.temRe;
176 double temIm[] = myfft.temIm;
177 //Permute() operation
178 permute(myfft, outputRe, outputIm, inputRe, inputIm);
180 //System.printString("ready to twiddle");
181 for (int factorIndex = 0; factorIndex < myfft.NumofFactors; factorIndex++)
182 twiddle(factorIndex, myfft, temRe, temIm, outputRe, outputIm);
184 //System.printString("ready to copy");
185 // Copy the output[] data to input[], so the output can be
186 // returned in the input array.
187 for (int i = 0; i < myfft.N; i++) {
188 inputRe[i] = outputRe[i];
189 inputIm[i] = outputIm[i];
193 private static void permute(fft1d myfft, double[] outputRe, double[] outputIm, double[] inputRe, double[] inputIm) {
194 int count[] = new int[myfft.MaxFactorsNumber];
198 for (int i = 0; i < myfft.N - 1; i++) {
199 outputRe[i] = inputRe[k];
200 outputIm[i] = inputIm[k];
202 k = k + myfft.remain[j];
203 count[0] = count[0] + 1;
204 while (count[j] >= myfft.factors[j]) {
210 tmp = myfft.remain[j - 1];
211 k = k - tmp + myfft.remain[j + 1];
213 count[j] = count[j] + 1;
216 outputRe[myfft.N - 1] = inputRe[myfft.N - 1];
217 outputIm[myfft.N - 1] = inputIm[myfft.N - 1];
218 } // End of function permute().
220 private static void twiddle(int factorIndex, fft1d myfft, double[] temRe, double[] temIm,
221 double[] outputRe, double[] outputIm) {
223 int sofarRadix = myfft.sofar[factorIndex];
224 int radix = myfft.factors[factorIndex];
225 int remainRadix = myfft.remain[factorIndex];
227 double tem; // Temporary variable to do data exchange.
229 double W = 2 * (double) Math.setPI() / (sofarRadix * radix);
230 double cosW = (double) Math.cos(W);
231 double sinW = -(double) Math.sin(W);
233 double twiddleRe[] = new double[radix];
234 double twiddleIm[] = new double[radix];
235 double twRe = 1.0f, twIm = 0f;
237 //Initialize twiddle addBk.address variables.
238 int dataOffset = 0, groupOffset = 0, address = 0;
240 for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
241 //System.printString("datano="+dataNo);
242 if (sofarRadix > 1) {
247 for (int i = 2; i < radix; i++) {
248 twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
249 twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
251 tem = cosW * twRe - sinW * twIm;
252 twIm = sinW * twRe + cosW * twIm;
255 for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
256 //System.printString("groupNo="+groupNo);
257 if ((sofarRadix > 1) && (dataNo > 0)) {
258 temRe[0] = outputRe[address];
259 temIm[0] = outputIm[address];
262 address = address + sofarRadix;
263 temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
264 twiddleIm[blockIndex] * outputIm[address];
265 temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
266 twiddleIm[blockIndex] * outputRe[address];
268 } while (blockIndex < radix);
270 for (int i = 0; i < radix; i++) {
271 //System.printString("temRe.length="+temRe.length);
272 //System.printString("i = "+i);
273 temRe[i] = outputRe[address];
274 temIm[i] = outputIm[address];
275 address += sofarRadix;
278 //System.printString("radix="+radix);
280 tem = temRe[0] + temRe[1];
281 temRe[1] = temRe[0] - temRe[1];
283 tem = temIm[0] + temIm[1];
284 temIm[1] = temIm[0] - temIm[1];
286 } else if( radix == 3) {
287 double t1Re = temRe[1] + temRe[2];
288 double t1Im = temIm[1] + temIm[2];
289 temRe[0] = temRe[0] + t1Re;
290 temIm[0] = temIm[0] + t1Im;
292 double m1Re = myfft.cos2to3PI * t1Re;
293 double m1Im = myfft.cos2to3PI * t1Im;
294 double m2Re = myfft.sin2to3PI * (temIm[1] - temIm[2]);
295 double m2Im = myfft.sin2to3PI * (temRe[2] - temRe[1]);
296 double s1Re = temRe[0] + m1Re;
297 double s1Im = temIm[0] + m1Im;
299 temRe[1] = s1Re + m2Re;
300 temIm[1] = s1Im + m2Im;
301 temRe[2] = s1Re - m2Re;
302 temIm[2] = s1Im - m2Im;
303 } else if(radix == 4) {
305 } else if(radix == 5) {
306 fft5(myfft, temRe, temIm);
307 } else if(radix == 8) {
308 fft8(myfft, temRe, temIm);
309 } else if(radix == 10) {
310 fft10(myfft, temRe, temIm);
312 fftPrime(radix, temRe, temIm);
314 address = groupOffset;
315 for (int i = 0; i < radix; i++) {
316 outputRe[address] = temRe[i];
317 outputIm[address] = temIm[i];
318 address += sofarRadix;
320 groupOffset += sofarRadix * radix;
321 address = groupOffset;
323 groupOffset = ++dataOffset;
324 address = groupOffset;
326 } //twiddle operation
328 // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
329 private static void fft4(double dataRe[], double dataIm[]) {
330 double t1Re,t1Im, t2Re,t2Im;
331 double m2Re,m2Im, m3Re,m3Im;
333 t1Re = dataRe[0] + dataRe[2];
334 t1Im = dataIm[0] + dataIm[2];
335 t2Re = dataRe[1] + dataRe[3];
336 t2Im = dataIm[1] + dataIm[3];
338 m2Re = dataRe[0] - dataRe[2];
339 m2Im = dataIm[0] - dataIm[2];
340 m3Re = dataIm[1] - dataIm[3];
341 m3Im = dataRe[3] - dataRe[1];
343 dataRe[0] = t1Re + t2Re;
344 dataIm[0] = t1Im + t2Im;
345 dataRe[2] = t1Re - t2Re;
346 dataIm[2] = t1Im - t2Im;
347 dataRe[1] = m2Re + m3Re;
348 dataIm[1] = m2Im + m3Im;
349 dataRe[3] = m2Re - m3Re;
350 dataIm[3] = m2Im - m3Im;
351 } // End of function fft4().
353 // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
354 private static void fft5(fft1d myfft, double dataRe[], double dataIm[]) {
355 double t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
356 double m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
357 double s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
359 t1Re = dataRe[1] + dataRe[4];
360 t1Im = dataIm[1] + dataIm[4];
361 t2Re = dataRe[2] + dataRe[3];
362 t2Im = dataIm[2] + dataIm[3];
363 t3Re = dataRe[1] - dataRe[4];
364 t3Im = dataIm[1] - dataIm[4];
365 t4Re = dataRe[3] - dataRe[2];
366 t4Im = dataIm[3] - dataIm[2];
370 dataRe[0] = dataRe[0] + t5Re;
371 dataIm[0] = dataIm[0] + t5Im;
373 m1Re = myfft.c51 * t5Re;
374 m1Im = myfft.c51 * t5Im;
375 m2Re = myfft.c52 * (t1Re - t2Re);
376 m2Im = myfft.c52 * (t1Im - t2Im);
377 m3Re = -(myfft.c53) * (t3Im + t4Im);
378 m3Im = myfft.c53 * (t3Re + t4Re);
379 m4Re = -(myfft.c54) * t4Im;
380 m4Im = myfft.c54 * t4Re;
381 m5Re = -(myfft.c55) * t3Im;
382 m5Im = myfft.c55 * t3Re;
388 s1Re = dataRe[0] + m1Re;
389 s1Im = dataIm[0] + m1Im;
395 dataRe[1] = s2Re + s3Re;
396 dataIm[1] = s2Im + s3Im;
397 dataRe[2] = s4Re + s5Re;
398 dataIm[2] = s4Im + s5Im;
399 dataRe[3] = s4Re - s5Re;
400 dataIm[3] = s4Im - s5Im;
401 dataRe[4] = s2Re - s3Re;
402 dataIm[4] = s2Im - s3Im;
403 } // End of function fft5().
405 private static void fft8(fft1d myfft, double[] temRe, double[] temIm) {
406 double data1Re[] = new double[4];
407 double data1Im[] = new double[4];
408 double data2Re[] = new double[4];
409 double data2Im[] = new double[4];
412 // To improve the speed, use direct assaignment instead for loop here.
413 data1Re[0] = temRe[0];
414 data2Re[0] = temRe[1];
415 data1Re[1] = temRe[2];
416 data2Re[1] = temRe[3];
417 data1Re[2] = temRe[4];
418 data2Re[2] = temRe[5];
419 data1Re[3] = temRe[6];
420 data2Re[3] = temRe[7];
422 data1Im[0] = temIm[0];
423 data2Im[0] = temIm[1];
424 data1Im[1] = temIm[2];
425 data2Im[1] = temIm[3];
426 data1Im[2] = temIm[4];
427 data2Im[2] = temIm[5];
428 data1Im[3] = temIm[6];
429 data2Im[3] = temIm[7];
431 fft4(data1Re, data1Im);
432 fft4(data2Re, data2Im);
434 tem = myfft.OnetoSqrt2 * (data2Re[1] + data2Im[1]);
435 data2Im[1] = myfft.OnetoSqrt2 * (data2Im[1] - data2Re[1]);
438 data2Im[2] = -data2Re[2];
440 tem = myfft.OnetoSqrt2 * (data2Im[3] - data2Re[3]);
441 data2Im[3] = -(myfft.OnetoSqrt2) * (data2Re[3] + data2Im[3]);
444 temRe[0] = data1Re[0] + data2Re[0];
445 temRe[4] = data1Re[0] - data2Re[0];
446 temRe[1] = data1Re[1] + data2Re[1];
447 temRe[5] = data1Re[1] - data2Re[1];
448 temRe[2] = data1Re[2] + data2Re[2];
449 temRe[6] = data1Re[2] - data2Re[2];
450 temRe[3] = data1Re[3] + data2Re[3];
451 temRe[7] = data1Re[3] - data2Re[3];
453 temIm[0] = data1Im[0] + data2Im[0];
454 temIm[4] = data1Im[0] - data2Im[0];
455 temIm[1] = data1Im[1] + data2Im[1];
456 temIm[5] = data1Im[1] - data2Im[1];
457 temIm[2] = data1Im[2] + data2Im[2];
458 temIm[6] = data1Im[2] - data2Im[2];
459 temIm[3] = data1Im[3] + data2Im[3];
460 temIm[7] = data1Im[3] - data2Im[3];
461 } // End of function fft8().
463 private static void fft10(fft1d myfft, double[] temRe, double[] temIm) {
464 double data1Re[] = new double[5];
465 double data1Im[] = new double[5];
466 double data2Re[] = new double[5];
467 double data2Im[] = new double[5];
469 // To improve the speed, use direct assaignment instead for loop here.
470 data1Re[0] = temRe[0];
471 data2Re[0] = temRe[5];
472 data1Re[1] = temRe[2];
473 data2Re[1] = temRe[7];
474 data1Re[2] = temRe[4];
475 data2Re[2] = temRe[9];
476 data1Re[3] = temRe[6];
477 data2Re[3] = temRe[1];
478 data1Re[4] = temRe[8];
479 data2Re[4] = temRe[3];
481 data1Im[0] = temIm[0];
482 data2Im[0] = temIm[5];
483 data1Im[1] = temIm[2];
484 data2Im[1] = temIm[7];
485 data1Im[2] = temIm[4];
486 data2Im[2] = temIm[9];
487 data1Im[3] = temIm[6];
488 data2Im[3] = temIm[1];
489 data1Im[4] = temIm[8];
490 data2Im[4] = temIm[3];
492 fft5(myfft, data1Re, data1Im);
493 fft5(myfft, data2Re, data2Im);
495 temRe[0] = data1Re[0] + data2Re[0];
496 temRe[5] = data1Re[0] - data2Re[0];
497 temRe[6] = data1Re[1] + data2Re[1];
498 temRe[1] = data1Re[1] - data2Re[1];
499 temRe[2] = data1Re[2] + data2Re[2];
500 temRe[7] = data1Re[2] - data2Re[2];
501 temRe[8] = data1Re[3] + data2Re[3];
502 temRe[3] = data1Re[3] - data2Re[3];
503 temRe[4] = data1Re[4] + data2Re[4];
504 temRe[9] = data1Re[4] - data2Re[4];
506 temIm[0] = data1Im[0] + data2Im[0];
507 temIm[5] = data1Im[0] - data2Im[0];
508 temIm[6] = data1Im[1] + data2Im[1];
509 temIm[1] = data1Im[1] - data2Im[1];
510 temIm[2] = data1Im[2] + data2Im[2];
511 temIm[7] = data1Im[2] - data2Im[2];
512 temIm[8] = data1Im[3] + data2Im[3];
513 temIm[3] = data1Im[3] - data2Im[3];
514 temIm[4] = data1Im[4] + data2Im[4];
515 temIm[9] = data1Im[4] - data2Im[4];
516 } // End of function fft10().
518 private static void fftPrime(int radix, double[] temRe, double[] temIm) {
520 double W = 2 * (double) Math.setPI() / radix;
521 double cosW = (double) Math.cos(W);
522 double sinW = -(double) Math.sin(W);
523 double WRe[] = new double[radix];
524 double WIm[] = new double[radix];
531 for (int i = 2; i < radix; i++) {
532 WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
533 WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
536 // FFT of prime length data, using DFT, can be improved in the future.
537 double rere, reim, imre, imim;
539 int max = (radix + 1) / 2;
541 double tem1Re[] = new double[max];
542 double tem1Im[] = new double[max];
543 double tem2Re[] = new double[max];
544 double tem2Im[] = new double[max];
546 for (j = 1; j < max; j++) {
547 tem1Re[j] = temRe[j] + temRe[radix - j];
548 tem1Im[j] = temIm[j] - temIm[radix - j];
549 tem2Re[j] = temRe[j] - temRe[radix - j];
550 tem2Im[j] = temIm[j] + temIm[radix - j];
553 for (j = 1; j < max; j++) {
556 temRe[radix - j] = temRe[0];
557 temIm[radix - j] = temIm[0];
559 for (int i = 1; i < max; i++) {
560 rere = WRe[k] * tem1Re[i];
561 imim = WIm[k] * tem1Im[i];
562 reim = WRe[k] * tem2Im[i];
563 imre = WIm[k] * tem2Re[i];
565 temRe[radix - j] += rere + imim;
566 temIm[radix - j] += reim - imre;
567 temRe[j] += rere - imim;
568 temIm[j] += reim + imre;
575 for (j = 1; j < max; j++) {
576 temRe[0] = temRe[0] + tem1Re[j];
577 temIm[0] = temIm[0] + tem2Im[j];
579 } // End of function fftPrime().