//l=0时，pr存放采样输入的实部，返回离散傅立叶变换的模；l=1时，pr存放傅立叶变换的实部，返回逆傅立叶变换的模
//l=0时，pi存放采样输入的虚部，返回离散傅立叶变换的幅角；l=1时，pi存放傅立叶变换的虚部，返回逆傅立叶变换的幅角
//l=0时，fr返回傅立叶变换的实部;l=1时返回逆傅立叶变换的实部
//l=0时，fi返回傅立叶变换的虚部;l=1时返回逆傅立叶变换的虚部 (l=0时正变换，l=1时逆变换)
//il=0时，表示不要求本函数计算傅立叶变换或逆傅立叶变换的模或幅角
//il=1时，表示要求本函数计算傅立叶变换或逆傅立叶变换的模或幅角
#include <math.h> //快速傅立叶变换子程序
void kfft(pr, pi, n, k, fr, fi, l, il)
int n, k, l, il;
float pr[], pi[], fr[], fi[];
{
int it, m, is, i, j, nv, l0;
float p, q, s, vr, vi, poddr, poddi;
for (it = 0; it < n; it++)
{
m = it;
is = 0;
for (i = 0; i < k; 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; 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; 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*nv;it = it+nv)
for (j = 0; j < nv/2; 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; i++)
{
fr[i] = fr[i]/(1.0*n);
fi[i] = fi[i]/(1.0*n);
}
if (il != 0)
for (i = 0; i < n; 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;
}
return;
}
