Path: utzoo!utgpu!jarvis.csri.toronto.edu!mailrus!tut.cis.ohio-state.edu!ucbvax!hplabs!hp-pcd!hplsla!jima
From: jima@hplsla.HP.COM (Jim Adcock)
Newsgroups: comp.lang.c++
Subject: Re: Vector and Matrix classes
Message-ID: <6590112@hplsla.HP.COM>
Date: 5 May 89 17:01:24 GMT
References: <1878@blake.acs.washington.edu>
Organization: HP Lake Stevens, WA
Lines: 60

#include <complex.h>
#include <stream.h>

// Fft a double precision complex array x of size (2 to the m)
// based on a fortran routine from Oppenheim and Schafer,
// though that's not the original source, I don't think.
// Even considering its simplicity, this fft does surprisingly well.

int 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;
        }
    }
    return m;
}

#define SZ 8
#define LGSZ 3

main()
{
    complex x[SZ];
    for (unsigned i = 0; i < SZ; ++i)
    {
	for (unsigned j = 0; j < SZ; ++j) x[j] = 0.0;
	x[i] = 1.0;
	fft(x, LGSZ);
	for (j = 0; j < SZ; ++j) cout << x[j] << " ";
	cout << "\n";
    }
}