/**********************fft programe*********************/
/*FFT-基2DIT-FFT算法*/
/*xr：=信号序列实部，xi：=信号序列虚部，N=2^M*/
/*计算结果X(k)的实部和虚部分别存储在数组xr和xi中*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
const float pp=3.14159;
void FFT(double xr[],double xi[],int N,int M)
{
int L,B,J,P,k,i; /*L级数，B=2^L-1为第L级B个不同旋转因子,P=j*2^M-L*/
float pi;
double ps;
double rx_kb,ix_kb; /*分别代表X(K+B)的实部和虚部*/
double rcf[64],icf[64]; /*rcf存储旋转因子实部，icf存储旋转因子虚部*/
/*旋转因子数组长根据需要调整*/
reverse_order(xr,xi,N) ; /*调用倒序子程序*/
/*计算各级蝶形*/
for(L=1;L<=M;L++)
{
pi=3.14159;
B=(int)(pow(2,(L-1))+0.5);
for(J=0;J<=B-1;J++) /*第L级B个不同的旋转因子*/
{
P=J*((int)(pow(2,(M-L))+0.5)); /*2^M-L个群体乘以Ｊ，共Ｎ/2个蝶形*/
ps=(2*pi/N)*P; /*计算旋转因子*/
rcf[P]=cos(ps);
icf[P]=-sin(ps);
for(k=J;k<=N-1;k+=(int)(pow(2,L)+0.5)) /*2^L为蝶距*/
{
rx_kb=xr[k+B]*rcf[P]-xi[k+B]*icf[P];
ix_kb=xi[k+B]*rcf[P]+xr[k+B]*icf[P];
xr[k+B]=xr[k]-rx_kb;
xi[k+B]=xi[k]-ix_kb;
xr[k]=xr[k]+rx_kb;
xi[k]=xi[k]+ix_kb;
}
}
}
}
/*倒序子程序*/
reverse_order(double xr[],double xi[],int N)
{
int LH,N1,I,J,K;
double T;
LH=N/2;
J=LH;
N1=N-2;
for(I=1;I<=N1;I++)
{
if(I<J)
{
T=xr[I]; xr[I]=xr[J]; xr[J]=T;
T=xi[I]; xi[I]=xi[J]; xi[J]=T;
}
K=LH;
while(J>=K)
{
J=J-K;
K=(int)(K/2+0.5);
}
J=J+K;
}
}
/**********************fft programe*********************/
/********************main programe**********************/
void main() /*以3cos(2pi*50t+pi/2)为例*/
{
int m,n;
double real[32],imag[32];
double result[32];
for(m=0;m<16;m++)
{
real[m]=3*cos(100*pp+pp/2+m*pp/8);
imag[m]=0;
}
FFT(real,imag,32,5);
clrscr();
for(m=0;m<16;m++)
{
/* printf("%8f",real[m]);
printf(" %f\n",imag[m]);*/
result[m]=sqrt(pow(real[m],2)+pow(imag[m],2));
printf("%8f\n",result[m]);
}
}
/********************main programe**********************/
