File Metadata

fft.c
View Options

	#include <math.h>

	#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

	void fft(double data, int nn, int ndim, int isign)
	{
	int i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
	int ibit, idim, k1, k2, n, nprev, nrem, ntot;
	double tempi, tempr;
	double theta, wi, wpi, wpr, wr, wtemp;

	ntot = 1;
	for (idim = 1; idim <= ndim; idim++)
	ntot *= nn[idim];
	nprev = 1;
	for (idim = ndim; idim >= 1; idim--)
	{
	n = nn[idim];
	nrem = ntot / (n * nprev);
	ip1 = nprev << 1;
	ip2 = ip1 * n;
	ip3 = ip2 * nrem;
	i2rev = 1;
	for (i2 = 1; i2 <= ip2; i2 += ip1)
	{
	if (i2 < i2rev)
	{
	for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2)
	{
	for (i3 = i1; i3 <= ip3; i3 += ip2)
	{
	i3rev = i2rev + i3 - i2;
	SWAP(data[i3], data[i3rev]);
	SWAP(data[i3+1], data[i3rev+1]);
	}
	}
	}
	ibit = ip2 >> 1;
	while (ibit >= ip1 && i2rev > ibit)
	{
	i2rev -= ibit;
	ibit >>= 1;
	}
	i2rev += ibit;
	}
	ifp1 = ip1;
	while (ifp1 < ip2)
	{
	ifp2 = ifp1 << 1;
	theta = isign * 6.28318530717959 / (ifp2 / ip1);
	wtemp = sin(0.5 * theta);
	wpr = -2.0 * wtemp * wtemp;
	wpi = sin(theta);
	wr = 1.0;
	wi = 0.0;
	for (i3 = 1; i3 <= ifp1; i3 += ip1)
	{
	for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2)
	{
	for (i2 = i1; i2 <= ip3; i2 += ifp2)
	{
	k1 = i2;
	k2 = k1 + ifp1;
	tempr = wr * data[k2] - wi * data[k2+1];
	tempi = wr * data[k2+1] + wi * data[k2];
	data[k2] = data[k1] - tempr;
	data[k2+1] = data[k1+1] - tempi;
	data[k1] += tempr;
	data[k1+1] += tempi;
	}
	}
	wr = (wtemp = wr) * wpr - wi * wpi + wr;
	wi = wi * wpr + wtemp * wpi + wi;
	}
	ifp1 = ifp2;
	}
	nprev *= n;
	}
	}

	#undef SWAP