首先说下为什么要FFT吧。主要进行的是时域和频域的变换。比如在设计滤波器的时候需要用到大量的卷积。这时候在时域上做卷积通过FFT转换到频域的时候，只需要对频域做乘加操作。再通过IFFT变换回时域即可。
#include "math.h"
#define PI
3.14159265358979
#define numlength
2048
void radix2(int n,short *xy,short *w);
void bitrev_cplx(int *x,short *index,int nx);
void digitrev_index(short *index, int n, int radix);
static short index1[64], w1[numlength];
static short Input[numlength*2];
static int Output[numlength];
static int radix;
int Num;
int i;
void main()
{
radix=2;
Num=numlength;
for (i=0;i
{
w1[2*i] = 32767*(-cos((double)i*2*PI/Num));
w1[2*i+1] = 32767*(-sin((double)i*2*PI/Num));
}
for (i=0;i
{
Input[2*i] = (short)((cos(PI*i/20.0)+cos(PI*i/10.0)+cos(PI*i/5.0))*0x80);
Input[2*i+1] = 0;
}
radix2(Num, Input, w1);
for (i=0; i
{
Output = sqrt(Input[2*i]*Input[2*i]+Input[2*i+1]*Input[2*i+1]);
}
digitrev_index(index1,Num,radix);
bitrev_cplx(Output, index1, Num);
}
//频率抽取FFT算法函数体
void radix2(int Num, short xy[], short w[])
{
short n1,n2,ie,ia,i,j,k,l;
short xt,yt,c,s;
n2 = Num;
ie = 1;
for (k=Num;k > 1;k = (k >> 1
) )
{
n1 = n2;
n2 = n2>>1;
ia = 0;
for (j=0; j < n2; j++)
{
c = w[2*ia];
s = w[2*ia+1];
ia = ia + ie;
for (i=j; i < Num; i += n1)
{
l = i + n2;
xt = xy[2*l] - xy[2*i];
xy[2*i] = xy[2*i] + xy[2*l];
yt = xy[2*l+1] - xy[2*i+1];
xy[2*i+1] = xy[2*i+1] + xy[2*l+1];
xy[2*l]
= (c*xt + s*yt)>>15;
xy[2*l+1] = (c*yt - s*xt)>>15;
}
}
ie = ie<<1;
}
}
//生成倒序表子程序
void digitrev_index(short *index, int n,int radix)
{
int
i, j, k;
short nbits, nbot, ntop, ndiff, n2, raddiv2;
nbits = 0;
i = n;
while (i > 1)
{
i = i >> 1;
nbits++;
}
raddiv2 = radix >> 1;
nbot = nbits >> raddiv2;
nbot = nbot << raddiv2-1;
ndiff
= nbits & raddiv2;
ntop
= nbot + ndiff;
n2
= 1 << ntop;
index[0] = 0;
for(i=1,j=n2/radix+1;i
{
index=j-1;
for (k = n2/radix; k*(radix-1) < j; k /= radix)
j -= k*(radix-1);
j += k;
}
index[n2-1] = n2 - 1;
}
//倒序子程序
void bitrev_cplx(int *x, short *index, int nx)
{
int
i;
short
i0, i1, i3;
short
j0, j1, j3;
int
xi0, xi1, xi3;
int
xj0, xj1, xj3;
short
t;
int
a, b, ia, ib, ibs;
int
mask;
int
nbits, nbot, ntop, ndiff, n2, halfn;
nbits = 0;
i = nx;
while (i > 1)
{
i = i >> 1;
nbits++;
}
nbot
= nbits >> 1;
ndiff
= nbits & 1;
ntop
= nbot + ndiff;
n2
= 1 << ntop;
mask
= n2 - 1;
halfn
= nx >> 1;
for (i0 = 0; i0 < halfn; i0 += 2)
{
b
= i0 & mask;
a
= i0 >> nbot;
if (!b) ia = index[a];
ib = index;
ibs = ib << nbot;
j0 = ibs + ia;
t
= i0 < j0;
xi0 = x[i0];
xj0 = x[j0];
if (t){x[i0] = xj0;
x[j0] = xi0;}
i1 = i0 + 1;
j1 = j0 + halfn;
xi1 = x[i1];
xj1 = x[j1];
x[i1] = xj1;
x[j1] = xi1;
i3 = i1 + halfn;
j3 = j1 + 1;
xi3 = x[i3];
xj3 = x[j3];
if (t){x[i3] = xj3;
x[j3] = xi3;}
}
}
| 欢迎光临 电子技术论坛_中国专业的电子工程师学习交流社区-中电网技术论坛 (http://bbs.eccn.com/) | Powered by Discuz! 7.0.0 |