Question about improving the efficiency of Fast Fourier Transform

[jeclow]

Hi, Can somebody please take a look at the following code and give me some hints on how to improve the efficiency of this code, like making it several times faster? I've tried to modify it many times, but none of them seemed to give me any significant improvement. Thank You very much!!

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

float pi=3.14159265358979;
#define matrix_element(i,j) ((i)*2+(j)-1)

int fft(int N,float *a,float *b)
{ int i;
float WNi[2];
float *x1,*x2,*X1,*X2;
if (N==1) {b[matrix_element(0,1)]=a[matrix_element(0,1)];b[matrix_element(0,2)]=a[matrix_element(0,2)];return 1;}

x1=(float *)malloc(N*sizeof(float));
x2=(float *)malloc(N*sizeof(float));
X1=(float *)malloc(N*sizeof(float));
X2=(float *)malloc(N*sizeof(float));
for(i=0;i<N;i++){
if (i%2){
x2[matrix_element((i/2),1)]=a[matrix_element(i,1)];
x2[matrix_element((i/2),2)]=a[matrix_element(i,2)];
}
else{
x1[matrix_element((i/2),1)]=a[matrix_element(i,1)];
x1[matrix_element((i/2),2)]=a[matrix_element(i,2)];
}
}

fft(N/2,x1,X1);
fft(N/2,x2,X2);
for(i=0;i<(N/2);i++){
WNi[matrix_element(0,1)]=cos(i*2*pi/N);
WNi[matrix_element(0,2)]=-sin(i*2*pi/N);
b[matrix_element(i,1)]=X1[matrix_element(i,1)]+X2[matrix_element(i,1)]*WNi[matrix_element(0,1)]-X2[matrix_element(i,2)]*WNi[matrix_element(0,2)];
b[matrix_element(i,2)]=X1[matrix_element(i,2)]+X2[matrix_element(i,1)]*WNi[matrix_element(0,2)]+X2[matrix_element(i,2)]*WNi[matrix_element(0,1)];
b[matrix_element((i+N/2),1)]=X1[matrix_element(i,1)]-X2[matrix_element(i,1)]*WNi[matrix_element(0,1)]+X2[matrix_element(i,2)]*WNi[matrix_element(0,2)];
b[matrix_element((i+N/2),2)]=X1[matrix_element(i,2)]-X2[matrix_element(i,1)]*WNi[matrix_element(0,2)]-X2[matrix_element(i,2)]*WNi[matrix_element(0,1)];
}
free(x1);
free(x2);
free(X1);
free(X2);
return 1;
}

int main()
{
FILE *fp0,*fp1;
clock_t start,finish;
float *input_array, *output_array;
float real,im;
int lengtha=1;
int i=0;
input_array=(float *)malloc(2*lengtha*sizeof(float));
fp0=fopen(&quot;infile&quot;,&quot;r&quot;
fscanf(fp0,&quot;%f&quot;,&real);
while(fscanf(fp0,&quot;%f&quot;,&im)!=EOF){
if((i+1)>lengtha){lengtha=lengtha*2;
if((input_array=(float *)realloc((void *)input_array,2*lengtha*sizeof(float)))==NULL) {printf(&quot;Memory allocation error\n&quot; return 0;}}
input_array[matrix_element(i,1)]=real;
input_array[matrix_element(i,2)]=im;
i++;
fscanf(fp0,&quot;%f&quot;,&real);
}
fclose(fp0);
output_array=(float *)malloc(2*lengtha*sizeof(float));
start=clock();
fft(lengtha,input_array,output_array);
finish=clock();
printf(&quot;Computed FFT in %f seconds&quot;,((float)(finish-start))/CLK_TCK);
fp1=fopen(&quot;outfile&quot;,&quot;w&quot;
for(i=0;i<=lengtha-1;i++){
fprintf(fp1,&quot;%f %f\n&quot;,output_array[matrix_element(i,1)],output_array[matrix_element(i,2)]);}
fclose(fp1);
}

 

[xwb]

It is basically the number of divisions and multiplications

  int allocsize = N*sizeof(float);
  x1=(float *)malloc(allocsize);
  x2=(float *)malloc(allocsize);
  X1=(float *)malloc(allocsize);  
  X2=(float *)malloc(allocsize);
  for(i=0;i<N;i++){
    int halfi = i << 1;
    if (i & 1){
      x2[matrix_element(halfi,1)]=a[matrix_element(i,1)];
      x2[matrix_element(halfi,2)]=a[matrix_element(i,2)];
    }
    else{
      x1[matrix_element(halfi,1)]=a[matrix_element(i,1)];
      x1[matrix_element(halfi,2)]=a[matrix_element(i,2)];
    }

But if you really want to speed it up significantly, use a lookup table for your sin and cos instead of calling the functions.

 

[xwb]

What does your data look like? I have tried some mods but I don't know how much faster it is and I don't know how big to make the sin/cos lookup tables.

 

[justice41]

You might also try looking for a non-recursive algorithm. For an average FT data set you are probably going to get into some pretty deep recursion. The danielson-lanczos algorithm presented in Numerical Recipes in C is non-recursive.

http://lib-www.lanl.gov/numerical/bookcpdf/c12-2.pdf

jaa

 

[jeclow]

XWB, justice41. Thank you so much for your reply.

I'm actaully expecting to improve it by 80%. If it can get the job done, the size of the program wouldn'y really matter. xwb, This is basically all my data. I'm trying to use your mod to see if there's big improvement or not. I'll let you know when I got the result. btw, do you think relocating memory would make the program faster?

Thanks again!!

 

[xwb]

It will be faster if you can use your own memory pool. malloc/free does some pretty strange things. As justice41 said, it would be faster if you can find a non-recursive algo.

 

[doherty]

on a side note check the bracketing and operator precedence
#define matrix_element(i,j) ((i)*2+(j)-1)
----------------------------------------
Friends are generally of the same sex, for when men and women agree,it is only in the conclusions; their reasons are always different.