Path: utzoo!attcan!uunet!lll-winken!ames!pasteur!ucbvax!hplabs!hp-pcd!hplsla!jima
From: jima@hplsla.HP.COM (Jim Adcock)
Newsgroups: comp.lang.c++
Subject: Re: Fast Fourier Transformer in C++
Message-ID: <6590086@hplsla.HP.COM>
Date: 23 Jan 89 23:25:17 GMT
References: <6250@thorin.cs.unc.edu>
Organization: HP Lake Stevens, WA
Lines: 46

/*******

For casual FFT'ers, here's a simple but elegant and efficient FFT
routine based on Rabiner and Gold pg 367

*******/

#include <complex.h>

// in-place fft a double precision complex array x of size (2 to the m)

void fft(complex* x, unsigned m)
{
    complex u, w, t;
    unsigned n = 1 << m;
    unsigned nv2 = n >> 1;
    unsigned nm1 = n - 1;
    unsigned j = 0;

    for (unsigned i = 0; i < nm1; ++i)
    {
        if (i < j) {t = x[j]; x[j] = x[i]; x[i] = t;}
        unsigned k = nv2;
        while ((k-1) < j) {j -= k; k >> 1;}
        j += k;
    }

    for (unsigned l = 1; l <= m; ++l)
    {
        unsigned le = 1 << l;
        unsigned le1 = le >> 1;
        complex u(1.0, 0.0);
        complex w(cos(M_PI/double(le1)),sin(M_PI/double(le1)));
        for (j = 0; j < le1; ++j)
        {
            for (i = j; i < n; i += le)
            {
                unsigned ip = i + le1;
                t = x[ip] * u;
                x[ip] = x[i] - t;
                x[i] += t;
            }
            u *= w;
        }
    }
}