96 lines
1.6 KiB
C
96 lines
1.6 KiB
C
/* FFT - Fast Fourier Transform */
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
#define FLOAT double
|
|
|
|
#define pi2 6.28318530717959
|
|
|
|
#define SWAP(a,b) { FLOAT temp = a; a = b; b = temp; }
|
|
|
|
static void four1 (FLOAT *data, int n)
|
|
{
|
|
int i, j, m, mmax;
|
|
|
|
i = 0;
|
|
j = 0;
|
|
while (i < n)
|
|
{
|
|
if (i < j)
|
|
{
|
|
SWAP(data[i], data[j]);
|
|
SWAP(data[i+1], data[j+1]);
|
|
}
|
|
m = n/2;
|
|
while (m >= 2 && j >= m)
|
|
{
|
|
j = j-m;
|
|
m = m/2;
|
|
}
|
|
i += 2;
|
|
j += m;
|
|
}
|
|
|
|
mmax = 2;
|
|
while (mmax < n)
|
|
{
|
|
FLOAT theta = pi2/mmax;
|
|
FLOAT x = sin (0.5*theta);
|
|
FLOAT wpr = -2.0*x*x;
|
|
FLOAT wpi= sin (theta);
|
|
FLOAT wr = 1.0;
|
|
FLOAT wi = 0.0;
|
|
FLOAT temp, tempr, tempi;
|
|
|
|
m = 0;
|
|
while (m < mmax)
|
|
{
|
|
i = m;
|
|
while (i < n)
|
|
{
|
|
j = i+mmax;
|
|
tempr = wr*data[j]-wi*data[j+1];
|
|
tempi = wi*data[j]+wr*data[j+1];
|
|
data[j] = data[i]-tempr;
|
|
data[j+1] = data[i+1]-tempi;
|
|
data[i] = data[i]+tempr;
|
|
data[i+1] = data[i+1]+tempi;
|
|
i = j+mmax;
|
|
}
|
|
temp = wr;
|
|
wr = temp*wpr-wi*wpi+wr;
|
|
wi = wi*wpr+temp*wpi+wi;
|
|
m += 2;
|
|
}
|
|
mmax *= 2;
|
|
}
|
|
}
|
|
|
|
#define N 1024
|
|
|
|
FLOAT data[N];
|
|
|
|
static FLOAT test ()
|
|
{
|
|
four1 (data, N);
|
|
return data[0];
|
|
}
|
|
|
|
int main (int argc, char *argv[])
|
|
{
|
|
int i;
|
|
FLOAT result;
|
|
|
|
for (i=0; i<N; i++)
|
|
data[i] = 0.0;
|
|
|
|
for (i=0; i<2000; i++)
|
|
result = test ();
|
|
|
|
if (result != 0.)
|
|
printf ("*** wrong result ***\n");
|
|
|
|
return 0;
|
|
}
|