www.pudn.com > fft_iccavr.rar > fft2.c

#include 
 
/*pr：输入实部 
pi：输入虚部 
n为2^k，n不小于数据个数 
l：为0表示正变换，1表示反变换 
il：1表示fr为取模结果。0表示fr输出实部fi输出虚部 
*/ 
void FFT(double* pr, double* pi, int n, int k, double* fr, double* fi, int l, int il) 
{ 
int it,m,is,i,j,nv,l0; 
double p,q,s,vr,vi,poddr,poddi; 
for(it = 0;it <= n - 1;it ++) 
{ 
m = it; 
is = 0; 
for(i = 0;i <= k - 1;i ++) 
{ 
j = m / 2; 
is = 2 * is + (m - 2 * j); 
m = j; 
} 
fr[it] = pr[is]; 
fi[it] = pi[is]; 
} 
pr[0] = 1.0; 
pi[0] = 0.0; 
p = 6.283185306 / (1.0 * n); 
pr[1] = cos(p); 
pi[1] = -sin(p); 
if(l != 0) 
pi[1] = -pi[1]; 
for(i = 2;i <= n - 1;i ++) 
{ 
p = pr[i - 1] * pr[1]; 
q = pi[i - 1] * pi[1]; 
s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]); 
pr[i] = p - q; 
pi[i] = s - p - q; 
} 
for(it = 0;it <= n - 2;it = it + 2) 
{ 
vr = fr[it]; 
vi = fi[it]; 
fr[it] = vr + fr[it + 1]; 
fi[it] = vi + fi[it + 1]; 
fr[it + 1] = vr - fr[it + 1]; 
fi[it + 1] = vi - fi[it + 1]; 
} 
m = n / 2; 
nv = 2; 
for(l0 = k - 2;l0 >= 0;l0 --) 
{ 
m = m / 2; 
nv = 2 * nv; 
for(it = 0;it <= (m - 1) * nv;it = it + nv) 
for(j = 0;j <= (nv / 2) - 1;j ++) 
{ 
p = pr[m * j] * fr[it + j + nv / 2]; 
q = pi[m * j] * fi[it + j + nv / 2]; 
s = pr[m * j] + pi[m * j]; 
s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]); 
poddr = p - q; 
poddi = s - p - q; 
fr[it + j + nv / 2] = fr[it + j] - poddr; 
fi[it + j + nv / 2] = fi[it + j] - poddi; 
fr[it + j] = fr[it + j] + poddr; 
fi[it + j] = fi[it + j] + poddi; 
} 
} 
if(l != 0) 
for(i = 0;i <= n - 1;i ++) 
{ 
fr[i] = fr[i] / (1.0 * n); 
fi[i] = fi[i] / (1.0 * n); 
} 
if(il != 0) 
for(i = 0;i <= n - 1;i ++) 
{ 
pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]); 
if(fabs(fr[i]) < 0.000001 * fabs(fi[i])) 
{ 
if((fi[i] * fr[i]) > 0) 
pi[i] = 90.0; 
else 
pi[i] = -90.0; 
} 
else 
pi[i] = atan(fi[i] / fr[i]) * 360.0 / 6.283185306; 
} 
}  
 
main() 
{ 
 while(1); 
}