Skip to content

Instantly share code, notes, and snippets.

@kristianlm kristianlm/fft.c
Created Dec 4, 2014

Embed
What would you like to do?
// 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);
}
@kristianlm

This comment has been minimized.

Copy link
Owner Author

kristianlm commented Dec 4, 2014

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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.