fft2.c

在ICCAVR 环境实现的fft,

#include<math.h>

/*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);
}