Zum Inhalt springen

Digitale bildgebende Verfahren: Transformationen/ Java-FFT

Aus Wikibooks

Klasse Complex

[Bearbeiten]

Die Java-Klasse Complex dient zur Speicherung von komplexen Zahlen mit zwei reellen Werten für den Realteil und für den Imaginärteil der komplexen Zahl sowie zur Berechnung des reellwertigen Betrags der komplexen Zahl:

/*
 * 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 or arguments
 */

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;
   }

   /**
   *  Computes the argument (phase angle with the positive real axis) of a complex number
   *  @return argument of complex number
   */
   public double argument ()
   {
      double cRe = this.re;
      double cIm = this.im;
      double result = java.lang.Math.atan2 (cRe, cIm);
      return result;
   }
}

Klasse FFT

[Bearbeiten]

Das folgende Programm in der Programmiersprache Java zur Berechnung von zweidimensionalen Fast-Fourier-Transformationen (FFT) kann mit dem Hauptprogramm FFT.main aufgerufen werden.

/*
 * 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];
			}
		}
	}

	/**
	 *  Example in one dimension
	 */
	private static void example1D ()
	{
		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");
	}

	/**
	 *  Example in two dimensions
	 */
	private static void example2D ()
	{
		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 ();
		}
	}

	/**
	 *  Main program with example calls for demonstration purposes
	 *  @param arguments: standard arguments
	 */
	public static void main (java.lang.String [] arguments)
	{
		example1D ();
		example2D ();
	}
}

/*
	Result table of fftArray after the call of example1D ()
	(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 example2D ()
	(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
*/

Ausgabe

[Bearbeiten]

Ausgabe der vom Programm berechneten Werte. Die gedoppelten Werte sind durchgestrichen und werden nach der Umsortierung (siehe Kommentar im Quelltext oben) nicht weiter verwendet:

8,4  0,4  8,0  0,4  0,4  0,4  0,4  0,4  0,4  0,4  0,4  0,4  0,4  0,4  8,0  0,4

31,8  2,0  0,2  0,2  0,2  0,2  0,2  2,0
 0,2  2,0  0,2  0,2  0,2  0,2  0,2  2,0
 1,4 16,2  1,4  1,4  1,4  1,4  1,4 16,2
 0,2  2,0  0,2  0,2  0,2  0,2  0,2  2,0
 0,2  2,0  0,2  0,2  0,2  0,2  0,2  2,0
 0,2  2,0  0,2  0,2  0,2  0,2  0,2  2,0
 1,4 16,2  1,4  1,4  1,4  1,4  1,4 16,2
 0,2  2,0  0,2  0,2  0,2  0,2  0,2  2,0