Small Ruby C extension to compute the FFT of an Array of power of 2 elements (Float).

FFT.c

/*
    fft extension to Array class.
    (C) Paolo Bosetti 2009 - University of Trento
    Credits to fft.c                                        
    (c) Douglas L. Jones                                  
    University of Illinois at Urbana-Champaign
    January 19, 1992
 */
#include "ruby.h"
#include <math.h>

/* this function is actually quite dirty and has to be clean up a lil */
void fft(int n, double x[], double y[])
{
    int m = (int)(log(n)/log(2));
    int i,j,k,n1,n2;
    double c,s,e,a,t1,t2;        

    j = 0; /* bit-reverse */
    n2 = n/2;
    for (i=1; i < n - 1; i++) {
        n1 = n2;
        while ( j >= n1 ) {
            j = j - n1;
            n1 = n1/2;
        }
        j = j + n1;

        if (i < j) {
            t1 = x[i];
            x[i] = x[j];
            x[j] = t1;
            t1 = y[i];
            y[i] = y[j];
            y[j] = t1;
        }
    }

    n1 = 0; /* FFT */
    n2 = 1;

    for (i=0; i < m; i++) {
        n1 = n2;
        n2 = n2 + n2;
        e = -6.283185307179586/n2;
        a = 0.0;

        for (j=0; j < n1; j++) {
            c = cos(a);
            s = sin(a);
            a = a + e;

            for (k=j; k < n; k=k+n2) {
                t1 = c*x[k+n1] - s*y[k+n1];
                t2 = s*x[k+n1] + c*y[k+n1];
                x[k+n1] = x[k] - t1;
                y[k+n1] = y[k] - t2;
                x[k] = x[k] + t1;
                y[k] = y[k] + t2;
            }
        }
    }
    return;
}

VALUE cArray;

static VALUE mFFT(VALUE self) {
    int n = (int)RARRAY(self)->len;
    float radix = log(n)/log(2);
    if (radix != (int)radix)
        rb_raise(rb_eRuntimeError, "Vector size must be power of 2 (was %d)", n);
    int i;
    double cAryRe[n], cAryIm[n];
    VALUE result = rb_ary_new2(n);
    VALUE complex[2];
    for (i=0; i<n; i++) {
        cAryRe[i] = NUM2DBL(rb_ary_entry(self, i));
        cAryIm[i] = 0.0;
    }
    fft(n,cAryRe,cAryIm);
    for (i=0; i<n; i++) {
        complex[0] = rb_float_new(cAryRe[i]);
        complex[1] = rb_float_new(cAryIm[i]);
        rb_ary_store(result, i, rb_class_new_instance(2,complex,rb_const_get(rb_cObject, rb_intern("Complex"))));
    }    
    return result;
}


void Init_FFT() {
    rb_require("complex");
    cArray = rb_define_class("Array", rb_cObject);
    rb_define_method(cArray, "fft", mFFT, 0);
}