Warning: Parameter 1 to Language::getMagic() expected to be a reference, value given in /home/wikija5/public_html/w/includes/StubObject.php on line 58
The Fast Fourier Transform in Java (part 1) - WikiJava
Saturday, 20th December 2014
Follow WikiJava on twitter now. @Wikijava

The Fast Fourier Transform in Java (part 1)

From WikiJava

Jump to: navigation, search
The author suggests:

buy this book


Here I'm going to show you a practical (and quite fast) Fast Fourier Transform I wrote in Java. Here I'll publish my comments on these codes and on the FFTs.

Contents

A foreword

"...the FFT has found application in almost every field of modern science, including fields as diverse as astronomy, acoustics, image processing, fluid dynamics, petroleum exploration, medicine and quantum physics, among many others. It is not an exercise in hyperbole to say that the world as we know it would be different without the FFT."

– review of David Bailey, to the book of Charles Van Loan "Computational Frameworks for the Fast Fourier Transform", March 1993.

The following code is based upon the Brigham's "The Fast Fourier Transform" (1975). Brigham describes in great detail the background you need to fully understand and write an optimized FFT. He even wrote an ALGOL60 and FORTRAN version of its algorithm. My humble implementation is heavily based on his great, and valuable, work, with a few modifications.

A quick and dirty explanation

There are a lot of different ways to compute this transformation. The slowest is the so-called DFT. In this particular FFT implementation we take 2 arrays of the same length (1 for the real part of the number, 1 for the imaginary part of the same number, so this implementation is referred as a "Complex one" - yep, there are some other kind of implementations, not complex) and we output 1 single array where every odd/even pair of elements are actually a complex number.

The Fast Fourier Transform (usually called FFT for short) is usually used to retrieve the "fundamental" frequency of a particular signal. With this transformation (and the relative magnitude) you could produce those beautiful graphs of the sound spectre.

The code

/**
 * @author Orlando Selenu
 * 
 */
public class FFTbase {
    /**
     * The Fast Fourier Transform (generic version, with NO optimizations).
     * 
     * @param inputReal
     *            an array of length n, the real part
     * @param inputImag
     *            an array of length n, the imaginary part
     * @param DIRECT
     *            TRUE = direct transform, FALSE = inverse transform
     * @return a new array of length 2n
     */
    public static double[] fft(final double[] inputReal, double[] inputImag,
	    boolean DIRECT) {
	// - n is the dimension of the problem
	// - nu is its logarithm in base e
	int n = inputReal.length;
 
	// If n is a power of 2, then ld is an integer (_without_ decimals)
	double ld = Math.log(n) / Math.log(2.0);
 
	// Here I check if n is a power of 2. If exist decimals in ld, I quit
	// from the function returning null.
	if (((int) ld) - ld != 0) {
	    System.out.println("The number of elements is not a power of 2.");
	    return null;
	}
 
	// Declaration and initialization of the variables
	// ld should be an integer, actually, so I don't lose any information in
	// the cast
	int nu = (int) ld;
	int n2 = n / 2;
	int nu1 = nu - 1;
	double[] xReal = new double[n];
	double[] xImag = new double[n];
	double tReal, tImag, p, arg, c, s;
 
	// Here I check if I'm going to do the direct transform or the inverse
	// transform.
	double constant;
	if (DIRECT)
	    constant = -2 * Math.PI;
	else
	    constant = 2 * Math.PI;
 
	// I don't want to overwrite the input arrays, so here I copy them. This
	// choice adds \Theta(2n) to the complexity.
	for (int i = 0; i < n; i++) {
	    xReal[i] = inputReal[i];
	    xImag[i] = inputImag[i];
	}
 
	// First phase - calculation
	int k = 0;
	for (int l = 1; l <= nu; l++) {
	    while (k < n) {
		for (int i = 1; i <= n2; i++) {
		    p = bitreverseReference(k >> nu1, nu);
		    // direct FFT or inverse FFT
		    arg = constant * p / n;
		    c = Math.cos(arg);
		    s = Math.sin(arg);
		    tReal = xReal[k + n2] * c + xImag[k + n2] * s;
		    tImag = xImag[k + n2] * c - xReal[k + n2] * s;
		    xReal[k + n2] = xReal[k] - tReal;
		    xImag[k + n2] = xImag[k] - tImag;
		    xReal[k] += tReal;
		    xImag[k] += tImag;
		    k++;
		}
		k += n2;
	    }
	    k = 0;
	    nu1--;
	    n2 /= 2;
	}
 
	// Second phase - recombination
	k = 0;
	int r;
	while (k < n) {
	    r = bitreverseReference(k, nu);
	    if (r > k) {
		tReal = xReal[k];
		tImag = xImag[k];
		xReal[k] = xReal[r];
		xImag[k] = xImag[r];
		xReal[r] = tReal;
		xImag[r] = tImag;
	    }
	    k++;
	}
 
	// Here I have to mix xReal and xImag to have an array (yes, it should
	// be possible to do this stuff in the earlier parts of the code, but
	// it's here to readibility).
	double[] newArray = new double[xReal.length * 2];
	double radice = 1 / Math.sqrt(n);
	for (int i = 0; i < newArray.length; i += 2) {
	    int i2 = i / 2;
	    // I used Stephen Wolfram's Mathematica as a reference so I'm going
	    // to normalize the output while I'm copying the elements.
	    newArray[i] = xReal[i2] * radice;
	    newArray[i + 1] = xImag[i2] * radice;
	}
	return newArray;
    }
 
    /**
     * The reference bitreverse function.
     */
    private static int bitreverseReference(int j, int nu) {
	int j2;
	int j1 = j;
	int k = 0;
	for (int i = 1; i <= nu; i++) {
	    j2 = j1 / 2;
	    k = 2 * k + j1 - 2 * j2;
	    j1 = j2;
	}
	return k;
    }
}

Comments from the users

To be notified via mail on the updates of this discussion you can login and click on watch at the top of the page


Comments on wikijava are disabled now, cause excessive spam.