C source code for Radix-2 FFT decimation-in-frequency algori

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define PI 3.14159265359
#define MAXPOW 24

struct complex
{
    double r;
    double i;
};

int pow_2[MAXPOW];
int pow_4[MAXPOW];

void twiddle(struct complex *W, int N, double stuff)
{
    W->r=cos(stuff*2.0*PI/(double)N);
    W->i=-sin(stuff*2.0*PI/(double)N);
}

void bit_reverse_reorder(struct complex *W, int N)
{
    int bits, i, j, k;
    double tempr, tempi;
    
    for (i=0; i<MAXPOW; i++)
	if (pow_2[i]==N) bits=i;

    for (i=0; i<N; i++)
    {
	j=0;
	for (k=0; k<bits; k++)
	    if (i&pow_2[k]) j+=pow_2[bits-k-1];

	if (j>i)  /** Only make "up" swaps */
	{
	    tempr=W[i].r;
	    tempi=W[i].i;
	    W[i].r=W[j].r;
	    W[i].i=W[j].i;
	    W[j].r=tempr;
	    W[j].i=tempi;
	}
    }
}
void bit_r4_reorder(struct complex *W, int N)
{
    int bits, i, j, k;
    double tempr, tempi;
    
    for (i=0; i<MAXPOW; i++)
	if (pow_2[i]==N) bits=i;

    for (i=0; i<N; i++)
    {
	j=0;
	for (k=0; k<bits; k+=2)
	{
	    if (i&pow_2[k]) j+=pow_2[bits-k-2];
	    if (i&pow_2[k+1]) j+=pow_2[bits-k-1];
	}

	if (j>i)  /** Only make "up" swaps */
	{
	    tempr=W[i].r;
	    tempi=W[i].i;
	    W[i].r=W[j].r;
	    W[i].i=W[j].i;
	    W[j].r=tempr;
	    W[j].i=tempi;
	}
    }
}

/** RADIX-4 FFT ALGORITHM */
void radix4(struct complex *x, int N)
{ 
    int    n2, k1, N1, N2;
    struct complex W, bfly[4];

    N1=4;
    N2=N/4;
    
    /** Do 4 Point DFT */ 
    for (n2=0; n2<N2; n2++)
    {
	/** Don't hurt the butterfly */
	bfly[0].r = (x[n2].r + x[N2 + n2].r + x[2*N2+n2].r + x[3*N2+n2].r);
	bfly[0].i = (x[n2].i + x[N2 + n2].i + x[2*N2+n2].i + x[3*N2+n2].i);

	bfly[1].r = (x[n2].r + x[N2 + n2].i - x[2*N2+n2].r - x[3*N2+n2].i);
	bfly[1].i = (x[n2].i - x[N2 + n2].r - x[2*N2+n2].i + x[3*N2+n2].r);

	bfly[2].r = (x[n2].r - x[N2 + n2].r + x[2*N2+n2].r - x[3*N2+n2].r);
	bfly[2].i = (x[n2].i - x[N2 + n2].i + x[2*N2+n2].i - x[3*N2+n2].i);

	bfly[3].r = (x[n2].r - x[N2 + n2].i - x[2*N2+n2].r + x[3*N2+n2].i);
	bfly[3].i = (x[n2].i + x[N2 + n2].r - x[2*N2+n2].i - x[3*N2+n2].r);


	/** In-place results */
	for (k1=0; k1<N1; k1++)
	{
	    twiddle(&W, N, (double)k1*(double)n2);
	    x[n2 + N2*k1].r = bfly[k1].r*W.r - bfly[k1].i*W.i;
	    x[n2 + N2*k1].i = bfly[k1].i*W.r + bfly[k1].r*W.i;
	}
    }
    
    /** Don't recurse if we're down to one butterfly */
    if (N2!=1)
	for (k1=0; k1<N1; k1++)
	{
	    radix4(&x[N2*k1], N2);
	}
}

/** RADIX-2 FFT ALGORITHM */
void radix2(struct complex *data, int N)
{
    int    n2, k1, N1, N2;
    struct complex W, bfly[2];

    N1=2;
    N2=N/2;
    
    /** Do 2 Point DFT */
    for (n2=0; n2<N2; n2++)
    {
	/** Don't hurt the butterfly */
	twiddle(&W, N, (double)n2);
	bfly[0].r = (data[n2].r + data[N2 + n2].r);
	bfly[0].i = (data[n2].i + data[N2 + n2].i);
	bfly[1].r = (data[n2].r - data[N2 + n2].r) * W.r - 
	    ((data[n2].i - data[N2 + n2].i) * W.i); 
	bfly[1].i = (data[n2].i - data[N2 + n2].i) * W.r +
	    ((data[n2].r - data[N2 + n2].r) * W.i);

	/** In-place results */
	for (k1=0; k1<N1; k1++)
	{
	    data[n2 + N2*k1].r = bfly[k1].r;
	    data[n2 + N2*k1].i = bfly[k1].i;
	}
    }
    
    /** Don't recurse if we're down to one butterfly */
    if (N2!=1)
	for (k1=0; k1<N1; k1++)
	    radix2(&data[N2*k1], N2);
}

void main(int argc, char *argv[])
{
    FILE   *infile;
    int    N, radix, numsamp;
    int    i;
    struct complex *data;
    double freq, phase, fs, A;
    int    dotime;
    struct timeval start, end;
    long   totaltime;
    
#ifdef GEN
    if (argc!=9)
    {
	printf("usage:\n");
	printf("    fft [A] [f] [phase] [fs] [num samp] [sequence length] [radix] [time]\n");
	printf("        output: DFT\n");
	exit(1);
    }
    
    
    sscanf(argv[1], "%lf", &A);
    sscanf(argv[2], "%lf", &freq);
    sscanf(argv[3], "%lf", &phase);
    sscanf(argv[4], "%lf", &fs);
    sscanf(argv[5], "%d", &numsamp);
    sscanf(argv[6], "%d", &N);
    sscanf(argv[7], "%d", &radix);
    sscanf(argv[8], "%d", &dotime);
#endif
#ifndef GEN
    if (argc<4)
    {
	printf("usage:\n");
	printf("    fft [input file] [sequence length] [radix]\n");
	printf("        output: DFT\n");
	exit(1);
    }
    else if ((infile=fopen(argv[1], "r"))==NULL)
    {
	printf("Error reading input sequence file: %s\n", argv[1]);
	exit(1);
    }
    
    sscanf(argv[2], "%d", &N);
    sscanf(argv[3], "%d", &radix);
    dotime=0;
#endif


    /** Set up power of two arrays */
    pow_2[0]=1;
    for (i=1; i<MAXPOW; i++)
	pow_2[i]=pow_2[i-1]*2;
    pow_4[0]=1;
    for (i=1; i<MAXPOW; i++)
	pow_4[i]=pow_4[i-1]*4;
    
    if ((data=malloc(sizeof(struct complex)*(size_t)N))==NULL)
    {
	fprintf(stderr, "Out of memory!\n");
	exit(1);
    }
    
    /** Generate cosine **/
#ifdef GEN
    for (i=0; i<N; i++)
    {
	data[i].r=0.0;
	data[i].i=0.0;
    }
    for (i=0; i<numsamp; i++)
	data[i].r=A*cos(2.0*PI*freq*i/fs - phase*PI/180);
#endif
#ifndef GEN
    for (i=0; i<N; i++)
    {
	fscanf(infile, "%lf", &data[i].r);
	data[i].i=0.0;
    }
#endif    

    gettimeofday(&start, (struct timezone *) 0);
    if (radix==2) radix2(data, N);
    if (radix==4) radix4(data, N);
    gettimeofday(&end, (struct timezone *) 0);
    totaltime=(end.tv_sec*1000000 + end.tv_usec)-(start.tv_sec*1000000
						  + start.tv_usec);
    if (radix==2) bit_reverse_reorder(data, N);
    if (radix==4) bit_r4_reorder(data, N);

    if (!dotime)
	for (i=0; i<N; i++)
	    printf("%f\n", sqrt(data[i].r*data[i].r +
				data[i].i*data[i].i));
    else
	printf("%ld us\n\n", totaltime);

#ifndef GEN
    fclose(infile);
#endif
}