fft.c

// code from http://paulbourke.net/miscellaneous/dft/
// help from http://www.codeproject.com/Articles/9388/How-to-implement-the-FFT-algorithm

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
/*
  This computes an in-place complex-to-complex FFT
  x and y are the real and imaginary arrays of 2^m points.
  dir =  1 gives forward transform
  dir = -1 gives reverse transform
*/

struct complex { double real, imag; } ;

short FFT(short int dir,long m, struct complex *buffer)
{
  long n,i,i1,j,k,i2,l,l1,l2;
  double c1,c2,tx,ty,t1,t2,u1,u2,z;

  /* Calculate the number of points */
  n = 1;
  for (i=0;i<m;i++)
    n *= 2;

  /* Do the bit reversal */
  i2 = n >> 1;
  j = 0;
  for (i=0;i<n-1;i++) {
    if (i < j) {
      tx = buffer[i].real;
      ty = buffer[i].imag;
      buffer[i].real = buffer[j].real;
      buffer[i].imag = buffer[j].imag;
      buffer[j].real = tx;
      buffer[j].imag = ty;
    }
    k = i2;
    while (k <= j) {
      j -= k;
      k >>= 1;
    }
    j += k;
  }

  /* Compute the FFT */
  c1 = -1.0;
  c2 = 0.0;
  l2 = 1;
  for (l=0;l<m;l++) {
    l1 = l2;
    l2 <<= 1;
    u1 = 1.0;
    u2 = 0.0;
    for (j=0;j<l1;j++) {
      for (i=j;i<n;i+=l2) {
        i1 = i + l1;
        t1 = u1 * buffer[i1].real - u2 * buffer[i1].imag;
        t2 = u1 * buffer[i1].imag + u2 * buffer[i1].real;
        buffer[i1].real = buffer[i].real - t1;
        buffer[i1].imag = buffer[i].imag - t2;
        buffer[i].real += t1;
        buffer[i].imag += t2;
      }
      z =  u1 * c1 - u2 * c2;
      u2 = u1 * c2 + u2 * c1;
      u1 = z;
    }
    c2 = sqrt((1.0 - c1) / 2.0);
    if (dir == 1)
      c2 = -c2;
    c1 = sqrt((1.0 + c1) / 2.0);
  }

  /* Scaling for forward transform */
  if (dir == 1) {
    for (i=0;i<n;i++) {
      buffer[i].real /= n;
      buffer[i].imag /= n;
    }
  }

  return(1);
}

// absolute value of complex number
double c_abs(struct complex n) {
  return n.real * n.real + n.imag * n.imag;
}

int main() {
  int m = 9;
  int i;
  int sample_rate = 8000; // default for arecord
  int len = pow(2, m);
  struct complex *buffer = malloc(sizeof(struct complex) * len);
  // measure in hertz!
  double factor = sample_rate / (double)len;

  for(i = 0 ; i < len ; i++) {
    double genHz = 1024;
    buffer[i].real = sin((i / (double)len) * genHz * M_PI);
    buffer[i].imag = 0;
  }

  fprintf(stderr, "len (bytes) = %d\nresolution (hz) = %.2f\n\n", len, factor);
  int quit = 0;
  while(!quit) {
    // fill buffer from stdin
    for(i = 0 ; i < len ; i++) {
      int ret = getchar();
      if(ret < 0) quit = 1;
      else buffer[i].real = ret;
      buffer[i].imag = 0;
    }
    // (/ 8000 2048.0)

    FFT(1, m, buffer);

    // find the fundamental frequency
    int max_index = -1;
    double max = 0;
    int freq_cutoff_low = 50. / factor; // ignore freqs below 50hz, not audible anyway

    for(i = 1 + freq_cutoff_low ; i < len/2 ; i++) {
      double v = c_abs(buffer[i]);
      if(v > max) {
        max_index = i;
        max = v;
      }
    }
    printf("peak freq (hz): %.0f \t(%f)\n", factor * max_index, max);
  }
  free(buffer);
}

I am playing a very annoying 2222hz tone from my speakers, picked up by my microphone. Run like this:

$ gcc -lm fft.c -o fft && arecord - | ./fft
len (bytes) = 512
resolution (hz) = 15.62

Recording WAVE '-' : Unsigned 8 bit, Rate 8000 Hz, Mono
peak freq (hz): 62  (37.748386)
peak freq (hz): 2219    (1.757652)
peak freq (hz): 2219    (1.741889)
peak freq (hz): 2219    (1.806563)
peak freq (hz): 2219    (1.925720)
peak freq (hz): 2219    (1.952391)
peak freq (hz): 2219    (1.991733)
peak freq (hz): 2219    (2.081786)