Robust/src/Benchmarks/Prefetch/2DFFT/dsm/fft1d.java

   1
   2 //Title:        1-d mixed radix FFT.
   3 //Version:
   4 //Copyright:    Copyright (c) 1998
   5 //Author:       Dongyan Wang
   6 //Company:      University of Wisconsin-Milwaukee.
   7 //Description:
   8 //  The number of DFT is factorized.
   9 //
  10 // Some short FFTs, such as length 2, 3, 4, 5, 8, 10, are used
  11 // to improve the speed.
  12 //
  13 // Prime factors are processed using DFT. In the future, we can
  14 // improve this part.
  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.
  19 //
  20 // A permute() function is used to make sure FFT can be calculated
  21 // in place.
  22 //
  23 // A triddle() function is used to perform the FFT.
  24 //
  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.
  29 //
  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.
  32 //
  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.
  37 //
  38
  39
  40 public class fft1d{
  41   // Maximum numbers of factors allowed.
  42   //private int MaxFactorsNumber = 30;
  43   public int MaxFactorsNumber;
  44
  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;
  50
  51   // TwotoFivePI   = 2*pi/5.
  52   // c51, c52, c53, c54, c55 are used in fft5().
  53   // c51 =(cos(TwotoFivePI)+cos(2*TwotoFivePI))/2-1.
  54   public double c51;
  55   // c52 =(cos(TwotoFivePI)-cos(2*TwotoFivePI))/2.
  56   public double c52;
  57   // c53 = -sin(TwotoFivePI).
  58   public double c53;
  59   // c54 =-(sin(TwotoFivePI)+sin(2*TwotoFivePI)).
  60   public double c54;
  61   // c55 =(sin(TwotoFivePI)-sin(2*TwotoFivePI)).
  62   public double c55;
  63
  64   // OnetoSqrt2 = 1/sqrt(2), used in fft8().
  65   public double OnetoSqrt2;
  66
  67   public int lastRadix;
  68
  69   int N;              // length of N point FFT.
  70   int NumofFactors;   // Number of factors of N.
  71   int maxFactor;      // Maximum factor of N.
  72
  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.
  76
  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;
  81
  82   // Constructor: FFT of Complex data.
  83   public fft1d(int N) {
  84     this.N = N;
  85     MaxFactorsNumber = 37;
  86     cos2to3PI = -1.5000f;
  87     sin2to3PI = 8.6602540378444E-01f;
  88     c51 = -1.25f;
  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;
  94     lastRadix = 0;
  95     maxFactor = 20;
  96     factorsWerePrinted = false;
  97     outputRe = global new double[N];
  98     outputIm = global new double[N];
  99
 100     factorize();
 101     //printFactors();
 102
 103     // Allocate memory for intermediate result of FFT.
 104     temRe = global new double[maxFactor]; //Check usage of this
 105     temIm = global new double[maxFactor];
 106   }
 107
 108   /*
 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");
 114       return;
 115     } else {
 116       this.inputRe = inputRe;
 117       this.inputIm = inputIm;
 118
 119       permute();
 120       //System.printString("ready to twiddle");
 121
 122       for (int factorIndex = 0; factorIndex < NumofFactors; factorIndex++)
 123         twiddle(factorIndex);
 124       //System.printString("ready to copy");
 125
 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];
 131       }
 132     }
 133   }
 134   */
 135
 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");
 142   }
 143
 144   public void factorize() {
 145     int radices[] = global new int[6];
 146     radices[0] = 2;
 147     radices[1] = 3;
 148     radices[2] = 4;
 149     radices[3] = 5;
 150     radices[4] = 8;
 151     radices[5] = 10;
 152     int temFactors[] = global new int[MaxFactorsNumber];
 153
 154     // 1 - point FFT, no need to factorize N.
 155     if (N == 1) {
 156       temFactors[0] = 1;
 157       NumofFactors = 1;
 158     }
 159
 160     // N - point FFT, N is needed to be factorized.
 161     int n = N;
 162     int index = 0;    // index of temFactors.
 163     int i = radices.length - 1;
 164
 165     while ((n > 1) && (i >= 0)) {
 166       if ((n % radices[i]) == 0) {
 167         n /= radices[i];
 168         temFactors[index++] = radices[i];
 169       } else
 170         i--;
 171     }
 172
 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)) {
 176       int test = 0;
 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.
 182           test = 1;
 183           //break;
 184         }
 185       }
 186     }
 187
 188     if (n > 1) {
 189       for (int k = 2; k < Math.sqrt(n) + 1; k++)
 190         while ((n % k) == 0) {
 191           n /= k;
 192           temFactors[index++] = k;
 193         }
 194       if (n > 1) {
 195         temFactors[index++] = n;
 196       }
 197     }
 198     NumofFactors = index;
 199     //if(temFactors[NumofFactors-1] > 10)
 200     //   maxFactor = n;
 201     //else
 202     //   maxFactor = 10;
 203
 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];
 208     }
 209
 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];
 216
 217     remain[0] = N / factors[0];
 218     sofar[0] = 1;
 219     for (i = 1; i < NumofFactors; i++) {
 220       sofar[i] = sofar[i - 1] * factors[i - 1];
 221       remain[i] = remain[i - 1] / factors[i];
 222     }
 223   }   // End of function factorize().
 224 /*
 225   private void permute() {
 226     int count[] = new int[MaxFactorsNumber];
 227     int j;
 228     int k = 0;
 229
 230     for (int i = 0; i < N - 1; i++) {
 231       outputRe[i] = inputRe[k];
 232       outputIm[i] = inputIm[k];
 233       j = 0;
 234       k = k + remain[j];
 235       count[0] = count[0] + 1;
 236       while (count[j] >= factors[j]) {
 237         count[j] = 0;
 238         k = k - (j == 0?N:remain[j - 1]) + remain[j + 1];
 239         j++;
 240         count[j] = count[j] + 1;
 241       }
 242     }
 243     outputRe[N - 1] = inputRe[N - 1];
 244     outputIm[N - 1] = inputIm[N - 1];
 245   }   // End of function permute().
 246   */
 247 /*
 248   private void twiddle(int factorIndex) {
 249     // Get factor data.
 250     int sofarRadix = sofar[factorIndex];
 251     int radix = factors[factorIndex];
 252     int remainRadix = remain[factorIndex];
 253
 254     double tem;   // Temporary variable to do data exchange.
 255
 256     double W = 2 * (double) Math.PI / (sofarRadix * radix);
 257     double cosW = (double) Math.cos(W);
 258     double sinW = -(double) Math.sin(W);
 259
 260     double twiddleRe[] = new double[radix];
 261     double twiddleIm[] = new double[radix];
 262     double twRe = 1.0f, twIm = 0f;
 263
 264     //Initialize twiddle addBk.address variables.
 265     int dataOffset = 0, groupOffset = 0, address = 0;
 266
 267     for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
 268       //System.printString("datano="+dataNo);
 269       if (sofarRadix > 1) {
 270         twiddleRe[0] = 1.0f;
 271         twiddleIm[0] = 0.0f;
 272         twiddleRe[1] = twRe;
 273         twiddleIm[1] = twIm;
 274         for (int i = 2; i < radix; i++) {
 275
 276
 277           twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
 278           twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
 279         }
 280         tem = cosW * twRe - sinW * twIm;
 281         twIm = sinW * twRe + cosW * twIm;
 282         twRe = tem;
 283       }
 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];
 289           int blockIndex = 1;
 290           do {
 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];
 296             blockIndex++;
 297           } while (blockIndex < radix);
 298         } else
 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;
 305           }
 306         //System.printString("radix="+radix);
 307         if(radix == 2) {
 308           case 2:
 309             tem = temRe[0] + temRe[1];
 310             temRe[1] = temRe[0] - temRe[1];
 311             temRe[0] = tem;
 312             tem = temIm[0] + temIm[1];
 313             temIm[1] = temIm[0] - temIm[1];
 314             temIm[0] = tem;
 315             break;
 316           case 3:
 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;
 321
 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;
 328
 329             temRe[1] = s1Re + m2Re;
 330             temIm[1] = s1Im + m2Im;
 331             temRe[2] = s1Re - m2Re;
 332             temIm[2] = s1Im - m2Im;
 333             break;
 334           case 4:
 335             fft4(temRe, temIm);
 336             break;
 337           case 5:
 338             fft5(temRe, temIm);
 339             break;
 340           case 8:
 341             fft8();
 342             break;
 343           case 10:
 344             fft10();
 345             break;
 346           default  :
 347             fftPrime(radix);
 348             break;
 349         }
 350         address = groupOffset;
 351         for (int i = 0; i < radix; i++) {
 352           outputRe[address] = temRe[i];
 353           outputIm[address] = temIm[i];
 354           address += sofarRadix;
 355         }
 356         groupOffset += sofarRadix * radix;
 357         address = groupOffset;
 358       }
 359       groupOffset = ++dataOffset;
 360       address = groupOffset;
 361     }
 362   } // End of function twiddle().
 363   */
 364 /*
 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;
 369
 370     t1Re = dataRe[0] + dataRe[2];
 371     t1Im = dataIm[0] + dataIm[2];
 372     t2Re = dataRe[1] + dataRe[3];
 373     t2Im = dataIm[1] + dataIm[3];
 374
 375     m2Re = dataRe[0] - dataRe[2];
 376     m2Im = dataIm[0] - dataIm[2];
 377     m3Re = dataIm[1] - dataIm[3];
 378     m3Im = dataRe[3] - dataRe[1];
 379
 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().
 389   */
 390 /*
 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;
 396
 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];
 405     t5Re = t1Re + t2Re;
 406     t5Im = t1Im + t2Im;
 407
 408     dataRe[0] = dataRe[0] + t5Re;
 409     dataIm[0] = dataIm[0] + t5Im;
 410
 411     m1Re = c51 * t5Re;
 412     m1Im = c51 * t5Im;
 413     m2Re = c52 * (t1Re - t2Re);
 414     m2Im = c52 * (t1Im - t2Im);
 415     m3Re = -c53 * (t3Im + t4Im);
 416     m3Im = c53 * (t3Re + t4Re);
 417     m4Re = -c54 * t4Im;
 418     m4Im = c54 * t4Re;
 419     m5Re = -c55 * t3Im;
 420     m5Im = c55 * t3Re;
 421
 422     s3Re = m3Re - m4Re;
 423     s3Im = m3Im - m4Im;
 424     s5Re = m3Re + m5Re;
 425     s5Im = m3Im + m5Im;
 426     s1Re = dataRe[0] + m1Re;
 427     s1Im = dataIm[0] + m1Im;
 428     s2Re = s1Re + m2Re;
 429     s2Im = s1Im + m2Im;
 430     s4Re = s1Re - m2Re;
 431     s4Im = s1Im - m2Im;
 432
 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().
 442   */
 443
 444   /*
 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];
 450     double tem;
 451
 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];
 461
 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];
 470
 471     fft4(data1Re, data1Im);
 472     fft4(data2Re, data2Im);
 473
 474     tem = OnetoSqrt2 * (data2Re[1] + data2Im[1]);
 475     data2Im[1] = OnetoSqrt2 * (data2Im[1] - data2Re[1]);
 476     data2Re[1] = tem;
 477     tem = data2Im[2];
 478     data2Im[2] = -data2Re[2];
 479     data2Re[2] = tem;
 480     tem = OnetoSqrt2 * (data2Im[3] - data2Re[3]);
 481     data2Im[3] = -OnetoSqrt2 * (data2Re[3] + data2Im[3]);
 482     data2Re[3] = tem;
 483
 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];
 492
 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().
 502   */
 503
 504   /*
 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];
 510
 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];
 522
 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];
 533
 534     fft5(data1Re, data1Im);
 535     fft5(data2Re, data2Im);
 536
 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];
 547
 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().
 559   */
 560
 561     /*
 562   public double sqrt(double d) {
 563     return Math.sqrt(d);
 564   }
 565   */
 566
 567     /*
 568   private void fftPrime(int radix) {
 569     // Initial WRe, WIm.
 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];
 575
 576     WRe[0] = 1;
 577     WIm[0] = 0;
 578     WRe[1] = cosW;
 579     WIm[1] = sinW;
 580
 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];
 584     }
 585
 586     // FFT of prime length data, using DFT, can be improved in the future.
 587     double rere, reim, imre, imim;
 588     int j, k;
 589     int max = (radix + 1) / 2;
 590
 591     double tem1Re[] = new double[max];
 592     double tem1Im[] = new double[max];
 593     double tem2Re[] = new double[max];
 594     double tem2Im[] = new double[max];
 595
 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];
 601     }
 602
 603     for (j = 1; j < max; j++) {
 604       temRe[j] = temRe[0];
 605       temIm[j] = temIm[0];
 606       temRe[radix - j] = temRe[0];
 607       temIm[radix - j] = temIm[0];
 608       k = j;
 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];
 614
 615         temRe[radix - j] += rere + imim;
 616         temIm[radix - j] += reim - imre;
 617         temRe[j] += rere - imim;
 618         temIm[j] += reim + imre;
 619
 620         k = k + j;
 621         if (k >= radix)
 622           k = k - radix;
 623       }
 624     }
 625     for (j = 1; j < max; j++) {
 626       temRe[0] = temRe[0] + tem1Re[j];
 627       temIm[0] = temIm[0] + tem2Im[j];
 628     }
 629   }   // End of function fftPrime().
 630   */
 631
 632 } // End of class FFT2d