Java implementation of the FFT algorithm

If you are looking use the FFT algorithm in a java application, you should probably take a look around for a well optimized and tested library. But if you want to get into it yourself, here is a start. I ported this code from C and I got that from “Numerical Recipes in C“. It was originally written by N. M. Brenner in C. This book is essential for anyone interested in DSP programming. I got it for a Numerical Methods class I took last year and I have still been finding uses for it.

For those of you who don’t know, the FFT algorithm allows you to view a discrete function (such as a sound wave) in the frequency domain. This means that you can take a buffer of sampled sound (or function) and retrieve the level of each frequency spectrum. This particular function returns a float array of the levels at each frequency bin b/w 0Hz and the Nyquist frequency. The first half of the array is the real part and the second is the imaginary. Here is a visual representation, the x axis is frequency bins, and the y axis is pressure (or energy):

This is not a very good picture, but, this is a 1 kHz sine wave. Here is the function:

public float[] four1(float data[], int nn, int isign) {

  int i, j, n, mmax, m, istep;
  float wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;

  n = nn << 1;
  j = 1;
  for(i = 1; i < n; i += 2) {
    if(j > i) {
      float temp;
      temp = data[j];
      data[j] = data[i];
      data[i] = temp;
      temp = data[j+1];
      data[j+1] = data[i+1];
      data[i+1] = temp;
    }
    m = n >> 1;
    while(m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
  mmax = 2;
  while(n > mmax) {
    istep = (mmax << 1);
    theta = isign*(6.28318530717959/mmax);
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;
    for(m = 1; m < mmax; m += 2) {
      for(i = m; i <= n; i += istep) {
        j = i+mmax;
        tempr = wr*data[j]-wi*data[j+1];
        tempi = wr*data[j+1]+wi*data[j];
        data[j] = data[i] - tempr;
        data[j+1] = data[i+1] - tempi;
        data[i] += tempr;
        data[i+1] += tempi;
      }
    wr = (wtemp=wr)*wpr-wi*wpi+wr;
    wi = wi*wpr+wtemp*wpi+wi;
    }
    mmax = istep;
  }
  return data;
}

data[] is the discrete function to be analyzed, nn is the number of complex points. This is generally data[].length/2, and isign is for inversion. Here is some code you can us to do this in Processing:

import krister.Ess.*;

AudioInput myInput;
float gain = 10.;
int bufferSize;

void setup() {
  size(512, 200);
  Ess.start(this);
  bufferSize=512;
  myInput = new AudioInput(bufferSize);
  myInput.gain(gain);
  frameRate(25);
  myInput.start();
  noFill();
  stroke(255, 255, 255);
  smooth();
}

void draw() {
  background(0, 0, 255);
  float[] buffer=four1(myInput.buffer2,
                       myInput.buffer2.length/4,1);
  stroke(255, 255, 255);
  for(int i=0; i
   line(i, height, i, height-(buffer[i]*1.4));
  }
  stroke(255, 0, 0);
  line(width/2, 0, width/2, height);
}

Of course, you need to drop the four1 function somewhere in there. This is a little crude for sound, but it will allow you to distinguish the difference in pitch of sounds with some accuracy.