Looking for a FFT in C/C++

[showtime]

I need FFT in c/c++ to deal with huge data series, 2^24 points. anyone has a high efficent fft in c or cpp?

 

[vadkudr]

Re: FFT in c/c++ needed

In DSPLib for TI C67xx there is C prototype sources for various type of FFT algorithms.

See ti.com

 

[fangf]

Re: FFT in c/c++ needed

Try fftw, www.fftw.org

 

[imperial_amazon]

Re: FFT in c/c++ needed

i recently did a program in C for N point FFT calculation using radix -2 DIF algoithm as a part of my lab course..

i'll post the program here.. please acknowledge

/* Author : Ragavan N B.E, D.E.C.E, */
/* Institute : Madras institute of Technology, Madras*/
/* Phone : 9840806628*/
/* e-mail : imperial_amazon@yahoo.com*/
/* In case you aren't able to follow the logic I used, feel free*/
/* to mail me.. i'll explain you in detail*/

# include <alloc.h>
# include <stdio.h>
# include <conio.h>
# include <math.h>

float pi=3.141592654;

struct Complex
{float real,imag;};

int BITREV(int N,int bininp) /*function to do the Bit reversal*/
{
int i,num,bitrevd=0;
num=log10(N)/log10(2);
for(i=0;i<(int)num;i++)
{if( ( (bininp<<i) & (N/2) )!=0) bitrevd+=(int)pow(2,i);}
return (bitrevd);
}


float wreal(int Npt,int i) /*calculate the real part re{exp(-j2*pi*n/N)}*/
{
if ( fabs(cos(2*pi*(float)i / (float)Npt)) <.0001 )
{return 0;}
else {return (float) ( cos(2*pi*(float)i / (float)Npt) ); }
}

float wimag(int Npt, int i) /*calculate the real part im{exp(-j2*pi*n/N)}*/
{
if ( fabs(-sin(2*pi*(float)i / (float)Npt)) < .0001)
{return 0;}
else
{return (float) ( -sin(2*pi*(float)i / (float)Npt) ); }
}



void decompose(int Npt, int N,int v,struct Complex (*x_time)[1024],struct Complex (*mainX)[1024])
{
struct Complex *x_tmp;int i;
x_tmp=malloc(2*Npt*sizeof(float));
for(i=0;i<Npt;i++) {x_tmp.real=0;x_tmp.imag=0;}
for(i=0;i< (N/2); i++)
{
auto float a,b,c,d;
x_tmp.real=x_time[v].real+x_time[v][i+(N/2)].real;
x_tmp.imag=x_time[v].imag+x_time[v][i+(N/2)].imag;
a=(x_time[v].real-x_time[v][i+(N/2)].real);
b=(x_time[v].imag-x_time[v][i+(N/2)].imag);
c=wreal(Npt,Npt*i/N);d=wimag(Npt,Npt*i/N);
x_tmp[i+(N/2)].real=( (a * c) - (b * d) );
x_tmp[i+(N/2)].imag=( (a * d) + (b * c) );
}
for(i=N*v;i<(N+(N*v));i++)
{
auto int s;
s=i-N*v;
if(x_tmp.real==0)
{mainX[0].real=abs(x_tmp.real);}
else
{mainX[0].real=x_tmp.real;}
if(x_tmp.imag==0)
{mainX[0].imag=abs(x_tmp.imag);}
else
{mainX[0].imag=x_tmp.imag;}
}
/* printf("\n\nDecomposed\n\n");
for(i=0;i<Npt;i++)
{
printf("\n%f \t %f",mainX[0].real,mainX[0].imag);
} */
free(x_tmp);
}

void main()
{
struct Complex (*x)[1024];
int i,Npt,N,stage,m,j,NumX,s;
clrscr();
printf("\n\nN = ");scanf("%d",&Npt);
stage=(int)(log10((double)Npt)/log10(2));
stage-=1;
x=malloc(2*2*Npt*sizeof(float));
for(i=0;i<Npt;i++)
{
x[0].real=0;
x[0].imag=0;
}
printf("\nx = ");
for(i=0;i<Npt;i++)
{scanf("%f", &x[0].real);}
for(i=0;i<Npt;i++)
{x[0].imag=0;}
decompose(Npt,Npt,0,x,x);
for(s=stage;s>0;s--)
{
struct Complex (*X1)[1024];
j=m=0;
X1=malloc(2*Npt*2*sizeof(float));
NumX=(int)(Npt/(int)pow(2,s));
for(i=0;i<NumX;i++)
{
for(;j<( m+((int)pow(2,s)) );j++)
{
(X1[j-m]).real=((float)((x[0][j]).real));
(X1[j-m]).imag=((float)((x[0][j]).imag));
}
m+=(int)pow(2,s);
}
for(i=0;i<NumX;i++)
{
decompose(Npt,(int)pow(2,s),i,X1,x);
}
free(X1);
}
/* printf("\n\nFFT output before Bit reversal\n\n");
for(i=0;i<Npt;i++)
{
printf("\n%f \t +j %f",x[0].real,x[0].imag);
} */
printf("\n\nResult of FFT\n\n");
for(i=0;i<Npt;i++)
{
printf("\n%f \t +j %f",x[0][BITREV(Npt,i)].real,x[0][BITREV(Npt,i)].imag);
}
free(x);
getch();

}

 

[Circuit_seller]

Re: FFT in c/c++ needed

Hi

For such a long FFT, radix-2 is not efficeint, try using radix-8 or more.

Regards

 

[sanbaba]

Re: FFT in c/c++ needed

// x is real coefficients
// y is imaginary coefficients

// By calling this function the resulting real and imaginary
// values will be stored in RealArray and ImaginaryArray

int main()
{

FFT( 1,Number1,RealArray,ImaginaryArray);// 1 for FFT
//Number1 Should be power of two
// Number1 = 5; for 32=2^5

FFT( -1,Number2,RealArray,ImaginaryArray);//-1 for IFFT

}

void FFT(int dir,int m,float *x,float *y)

{

int n,i,i1,j,k,i2,l,l1,l2;

float c1,c2,tx,ty,t1,t2,u1,u2,z;

//No of Pts



n = 1;

for (i=0;i<m;i++)

n *= 2;
/*

0000 -----> 0000 0

0001 -----> 1000 8

0010 -----> 0100 4

0011 -----> 1100 12

0100 -----> 0010 2

0101 -----> 1010

0110 -----> 0110

0111 -----> 1110

1000 -----> 0001

1001 -----> 1001

1010 -----> 0101

1011 -----> 1101

1100 -----> 0011

1101 -----> 1011

1110 -----> 0111

1111 -----> 1111

*/
i2 = n >> 1; // 2^n

j = 0;

for (i=0;i<n-1;i++)

{

if (i < j)

{

tx = x;

ty = y;

x = x[j];

y = y[j];

x[j] = tx;

y[j] = ty;

}


k = i2;

while (k <= j)

{

j -= k;

k >>= 1;

}

j += k;
}


c1 = -1;

c2 = 0;

l2 = 1;


for (l=0;l<m;l++)

{

l1 = l2;

l2 <<= 1;

u1 = 1;

u2 = 0;


for (j=0;j<l1;j++)

{

for (i=j;i<n;i+=l2)

{

i1 = i + l1;

t1 = (u1 * x[i1] - u2 * y[i1]);

t2 = (u1 * y[i1] + u2 * x[i1]);

x[i1] = x - t1;

y[i1] = y - t2;

x += t1;

y += t2;

}

z = ( u1 * c1 - u2 * c2);

u2 =(u1 * c2 + u2 * c1);

u1 = z;

}

c2 = sqrt((1.0 - c1) / 2.0);

if (dir == 1)

c2 = -c2;

c1 = sqrt((1.0 + c1) / 2.0);

}





if (dir == -1)

{

for (i=0;i<n;i++)

{

x = x/n;


y = y/n;


}

}

}