Fast Fourier Transform in C

The following code is useful to compute Frequency Spectrum from time series input. The input is an ascii file input.txt and the output is automatically displayed by GNUPLOT. You may need to modify sr (sampling rate in time domain (in seconds)) in the following code. The original code is from Numerical Recipes and some sources available in internet. 

The command is:
gcc -o fft fft.c -lm

./fft input.txt


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define PI    M_PI   
#define TWOPI    (2.0*PI)
void four1(float data[], int nn, int isign)
{
    int n, mmax, m, j, istep, i;
    float wtemp, wr, wpr, wpi, wi, theta;
    float tempr, tempi;
    n = nn << 1;  
    j = 1;
    for (i = 1; i < n; i += 2) {
    if (j > i) {
    tempr = data[j];       
    data[j] = data[i];   
    data[i] = tempr;       
    tempr = data[j+1];     
    data[j+1] = data[i+1]; 
    data[i+1] = tempr;    
    }
    m = n >> 1;
    while (m >= 2 && j > m) {
        j -= m;
        m >>= 1; 
    }
    j += m;  
    }
    mmax = 2;
    while (n > mmax) {
    istep = 2*mmax;
    theta = TWOPI/(isign*mmax); 
    wtemp = sin(0.5*theta); 
    wpr = -2.0*wtemp*wtemp; 
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;
    for (m = 1; m < mmax; m += 2) {
        for (i = m; i <= n; i += istep) {
        j =i + mmax;
        tempr = wr*data[j]   - wi*data[j+1];
        tempi = wr*data[j+1] + wi*data[j];
        data[j]   = data[i]   - tempr;
        data[j+1] = data[i+1] - tempi;
        data[i] += tempr;
        data[i+1] += tempi;
        }
        wr = (wtemp = wr)*wpr - wi*wpi + wr;
        wi = wi*wpr + wtemp*wpi + wi;

    }
    mmax = istep;
    }
}


int main(int argc, char * argv[])
{

int i;
int Nx;
int NFFT;
float *x;
float *X;
float f, f1, sr;
sr =0.001; //sampling rate in time domain

//count how many rows in input
FILE* fileku = fopen(argv[argc-1],"r");
int ch, Mx = 0;
do
{
ch = fgetc(fileku);
if(ch == '\n')
Mx++;
} while (ch != EOF);
if(ch != '\n' && Mx != 0)
Mx++;
fclose(fileku);
Nx=Mx-1;  // no of rows

//read input file
FILE* myfile = fopen(argv[argc-1],"r");
x = (float *) malloc(Nx * sizeof(float));
for(i=0;i<Nx;++i)
{
fscanf(myfile,"%f",&x[i]);
}
fclose(myfile);
NFFT = (int)pow(2.0, ceil(log((float)Nx)/log(2.0)));
/* allocate memory for NFFT complex numbers (note the +1) */
X = (float *) malloc((2*NFFT+1) * sizeof(float));
for(i=0; i<Nx; i++)
{
X[2*i+1] = x[i];
X[2*i+2] = 0.0;
}
for(i=Nx; i<NFFT; i++)
{
X[2*i+1] = 0.0;
X[2*i+2] = 0.0;
}
/* calculate FFT */
four1(X, NFFT, 1);



/* PLot result with GNUPLOT*/
    char * commandsForGnuplot[] = {"set title \"FFT RESULT\"", "set xlabel 'Frequency (Hz)'","set ylabel 'Amplitude Relative'","plot 'result' using 1:2 with lines"};
    FILE * temp = fopen("result", "w");
    FILE * gnuplotPipe = popen ("gnuplot -persistent", "w");
    f = 1;
    for(i=0; i<(NFFT/2)+1; i++)
    {
    f = i;
    f1=(1/(2*sr))*(f/(NFFT/2));
    fprintf(temp, "%2f  %.2f \n", f1,sqrt((X[2*i+1]*X[2*i+1]) + (X[2*i+2]*X[2*i+2])));
    }
    for (i=0; i < 4; i++)
    {
    fprintf(gnuplotPipe, "%s \n", commandsForGnuplot[i]);
    }
    return 0;
}