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 = new double[N];
  98     outputIm = new double[N];
  99
 100     factorize();
 101     //printFactors();
 102
 103     // Allocate memory for intermediate result of FFT.
 104     temRe = new double[maxFactor]; //Check usage of this
 105     temIm = new double[maxFactor];
 106   }
 107
 108   public void printFactors() {
 109     if (factorsWerePrinted) return;
 110     factorsWerePrinted = true;
 111     System.printString("factors.length = " + factors.length + "\n");
 112     for (int i = 0; i < factors.length; i++)
 113       System.printString("factors[i] = " + factors[i] + "\n");
 114   }
 115
 116   public void factorize() {
 117     int radices[] = new int[6];
 118     radices[0] = 2;
 119     radices[1] = 3;
 120     radices[2] = 4;
 121     radices[3] = 5;
 122     radices[4] = 8;
 123     radices[5] = 10;
 124     int temFactors[] = new int[MaxFactorsNumber];
 125
 126     // 1 - point FFT, no need to factorize N.
 127     if (N == 1) {
 128       temFactors[0] = 1;
 129       NumofFactors = 1;
 130     }
 131
 132     // N - point FFT, N is needed to be factorized.
 133     int n = N;
 134     int index = 0;    // index of temFactors.
 135     int i = radices.length - 1;
 136
 137     while ((n > 1) && (i >= 0)) {
 138       if ((n % radices[i]) == 0) {
 139         n /= radices[i];
 140         temFactors[index++] = radices[i];
 141       } else
 142         i--;
 143     }
 144
 145     // Substitute 2x8 with 4x4.
 146     // index>0, in the case only one prime factor, such as N=263.
 147     if ((index > 0) && (temFactors[index - 1] == 2)) {
 148       int test = 0;
 149       for (i = index - 2; (i >= 0) && (test == 0); i--) {
 150         if (temFactors[i] == 8) {
 151           temFactors[index - 1] = temFactors[i] = 4;
 152           // break out of for loop, because only one '2' will exist in
 153           // temFactors, so only one substitutation is needed.
 154           test = 1;
 155         }
 156       }
 157     }
 158
 159     if (n > 1) {
 160       for (int k = 2; k < Math.sqrt(n) + 1; k++)
 161         while ((n % k) == 0) {
 162           n /= k;
 163           temFactors[index++] = k;
 164         }
 165       if (n > 1) {
 166         temFactors[index++] = n;
 167       }
 168     }
 169     NumofFactors = index;
 170
 171     // Inverse temFactors and store factors into factors[].
 172     factors = new int[NumofFactors];
 173     for (i = 0; i < NumofFactors; i++) {
 174       factors[i] = temFactors[NumofFactors - i - 1];
 175     }
 176
 177     // Calculate sofar[], remain[].
 178     // sofar[]  : finished factors before the current stage.
 179     // factors[]: factors of N processed in the current stage.
 180     // remain[] : finished factors after the current stage.
 181
 182     sofar = new int[NumofFactors];
 183     remain = new int[NumofFactors];
 184
 185     remain[0] = N / factors[0];
 186     sofar[0] = 1;
 187     for (i = 1; i < NumofFactors; i++) {
 188       sofar[i] = sofar[i - 1] * factors[i - 1];
 189       remain[i] = remain[i - 1] / factors[i];
 190     }
 191   }   // End of function factorize().
 192 } // End of class FFT1d