Digitale bildgebende Verfahren: Transformationen/ Java-FFT

/*
 * Title: Complex
 * Last edit: 24 October 2025
 * Author: Markus Bautsch, Berlin
 * Programming Language: Java
 */

/**
 * This program is managing complex numbers and computes their absolute values
 */

public class Complex
{
   /**
    * real part of complex number
    */
   public double re;
   /**
    * imaginary part of complex number
    */
   public double im;

   /**
    * Constructor for the instance of a complex number
    * @param re: real part of complex number
    * @param im: imaginary part of complex number
    */
   public Complex (double re, double im)
   {
      this.re = re;
      this.im = im;
   }

   /**
    *  Computes the absolute value of a complex number
    *  @return absolute value of complex number
    */
   public double amplitude ()
   {
      double cRe = this.re;
      double cIm = this.im;
      double result = java.lang.Math.sqrt (cRe * cRe + cIm * cIm);
      return result;
   }
}

/*
 * Title: Fast-Fourier-Transformation, FFT
 * Last edit: 24 October 2025
 * Author: Markus Bautsch, Berlin
 * Programming Language: Java, translated from Component Pascal (4 January 2012)
 * Reference: Elbert Oran, FFT-Anwendungen, Oldenbourg, München, Wien, 1997
 */

/**
 * This program is carrying out fast Fourier transforms on one or two-dimensional complex arrays
 */

public class FFT
{
   /**
    *  Class variables
    */
   private int n;  // size of FFT array
   private int nu; // power of two of size of arrays
   private double [] cosArray;
   private double [] sinArray;
   private int    [] bitrev;

   /**
    *  Constructor creates and initialises arrays for computing FFT
    *  @param n: size of complex array
    */
   public FFT (int n)
   {
      this.nu = java.lang.Integer.numberOfTrailingZeros (n); // integer exponent of power of two of n
      if (this.nu > 1)
      {
         this.n = n;
         this.cosArray = new double [n];
         this.sinArray = new double [n];
         this.bitrev = new int [n];
         this.initArrays ();
      }
   }

   /**
    *  Initialises arrays for computing FFT
    */
   private void initArrays ()
   {
      double pi2N = 2 * java.lang.Math.PI / this.n;
      int i = this.n;
      do
      {
         i--;
         int bitset = 0; // empty set
         int nu1 = 0;
         int j = this.nu;
         do
         {
            j--;
            if ((i & (1 << j)) != 0) // ask for bit j in i
            {
               bitset = bitset | (1 << nu1); // set bit nu1 in bitset
            }
            nu1++;
         } while (j > 0);
         this.bitrev [i] = bitset;

         double arg = pi2N * i;
         this.cosArray [i] = java.lang.Math.cos (arg);
         this.sinArray [i] = java.lang.Math.sin (arg);
      } while (i > 0);
   }

   /**
    *  Sorts one-dimensional complex array according to FFT bit reversal algorithm
    *  @param h: one-dimensional complex array
    */
   private void sort (Complex [] h)
   {
      int l = h.length - 1;
      int k = 0;
      do
      {
         k++;
         int i = this.bitrev [k];
         if (i > k)
         {
            Complex temp = h [k];
            h [k] = h [i];
            h [i] = temp;
         }
      } while (k < l);
   }

   /**
    *  Calculates one-dimensional discrete Fourier transformation H of a complex array h.
    *  The result is stored in the same array h because of efficiency reasons.
    *  The two halfs of the result are not exchanged, i.e. the constant level of the Fourier transform is contained in h [0] !
    *  @param h: one-dimensional complex array
    */
   private void fft (Complex [] h)
   {
      h [0].re = 0.5 * (h [0].re + h [n - 1].re);
      h [0].im = 0.5 * (h [0].im + h [n - 1].im);
      int n2 = this.n / 2;
      for (int l = this.nu - 1; l >= 0; l--)
      {
         int bitset = (1 << l); // assign bit l;
         int k = 0;
         while (k < this.n) 
         {
            int kn2 = k + n2;
            for (int i = 0; i < n2; i++)
            {
               int p = this.bitrev [k / bitset];
               double cos = this.cosArray [p];
               double sin = this.sinArray [p];
               Complex hc = h [kn2];
               double hcre = hc.re;
               double hcim = hc.im;
               double tempRe = hcre * cos + hcim * sin;
               double tempIm = hcim * cos - hcre * sin;

               hc = h [k];
               hcre = hc.re;
               hcim = hc.im;
               h [kn2].re = hcre - tempRe;
               h [kn2].im = hcim - tempIm;
               h [k].re   = hcre + tempRe;
               h [k].im   = hcim + tempIm;
               k++;
               kn2++;
            }
            k = k + n2;
         }
         n2 = n2 / 2;
      }
      this.sort (h);
   }

   /**
    *  Calculates two-dimensional discrete Fourier transformation H2 of a complex array h2.
    *  The result is stored in the same array h2 because of efficiency reasons.
    *  @param h2: two-dimensional complex array
    */
   private void fft2D (Complex [] [] h2)
   {
      /* calc columns */
      for (int j = 0; j < this.n; j++)
      {
         this.fft (h2 [j]);
      }

      /* calc rows: */
      Complex [] ht = new Complex [this.n];
      for (int i = 0; i < this.n; i++)
      {
         for (int j = 0; j < this.n; j++)
         {
            ht [j] = h2 [j] [i];
         }
         this.fft (ht); // transform row
         for (int j = 0; j < this.n; j++)
         {
            h2 [j] [i] = ht [j];
         }
      }
   }

   /**
    *  Main program with example calls for demonstration purposes
    *  @param arguments: standard arguments
    */
   public static void main (java.lang.String [] arguments)
   {
      // one-dimensional example
      int n1d = 16;
      FFT fft1D = new FFT (n1d);
      Complex [] fftArray1D = new Complex [fft1D.n] ; // 1D-array
      for (int i = 0; i < fft1D.n; i++)
      {
         double re =
               java.lang.Math.sin (4 * java.lang.Math.PI * i / fft1D.n) - 0.5;
         double im = 0;
         fftArray1D [i] = new Complex (re, im);
      }
      fft1D.fft (fftArray1D);
      for (int i = 0; i < fft1D.n; i++)
      {
         java.lang.System.out.printf ("%4.1f", fftArray1D [i].amplitude ());
         java.lang.System.out.print (" ");
      }
      java.lang.System.out.println ("\n");

      // two-dimensional example
      int n2d = 8;
      FFT fft2D = new FFT (n2d);
      Complex [] [] fftArray2D = new Complex [fft2D.n] [fft2D.n]; // 2D-array must be square
      for (int j = 0; j < fft2D.n; j++)
      {
         for (int i = 0; i < fft2D.n; i++)
         {
            double re =
                  java.lang.Math.sin (2 * java.lang.Math.PI * i / fft2D.n) *
                  java.lang.Math.sin (4 * java.lang.Math.PI * j / fft2D.n) - 0.5;
            double im = 0;
            fftArray2D [i] [j] = new Complex (re, im);
         }
      }
      fft2D.fft2D (fftArray2D);
      for (int j = 0; j < fft2D.n; j++)
      {
         for (int i = 0; i < fft2D.n; i++)
         {
            java.lang.System.out.printf ("%4.1f", fftArray2D [i] [j].amplitude ());
            java.lang.System.out.print (" ");
         }
         java.lang.System.out.println ();
      }
   }
}

/*
   Result table of fftArray after the call of main
   (the two half sub blocks have to be exchanged and first number has been removed):
   0.4  0.4  0.4  0.4  0.4  8.0  0.4  8.4  0.4  8.0  0.4  0.4  0.4  0.4  0.4
 
   Result table of fftArray2D after the call of main
   (the four quarter sub blocks have to be exchanged diagonally and first row 0 as well as first column 0 have been removed):
   0.2 0.2  2.0  0.2  2.0 0.2 0.2
   1.4 1.4 16.2  1.4 16.2 1.4 1.4
   0.2 0.2  2.0  0.2  2.0 0.2 0.2
   0.2 0.2  2.0 31.8  2.0 0.2 0.2
   0.2 0.2  2.0  0.2  2.0 0.2 0.2
   1.4 1.4 16.2  1.4 16.2 1.4 1.4
   0.2 0.2  2.0  0.2  2.0 0.2 0.2
*/