fft.c

FFT

#define PI 3.14159265358979
#define M 10
#define N 1024
#include<math.h>
#include<stdio.h>
void radix2(int*xr,int*xi,short*wr,short*wi);
void ChangeOrder(int*xr,int*xi);
int xr[N],xi[N];
int x[2*N];
int xm[N];
int xlr[N],xli[N];
short w1[N],w2[N];
int y[2*N];
int ym[N];
int radix;
int nxl;
int i;
main()
{
double delta;
radix=2;
nxl=N;
delta=2*PI/nxl;
for(i=0;i<nxl/2;i++)
{
w1[i]=32767*(cos(i*delta));
w2[i]=32767*(sin(i*delta));
}
for(i=0;i<nxl;i++)
{
xr[i]=xlr[i]=(int)((cos(PI*i/5.0))*0x80);
xi[i]=xli[i]=1;
x[2*i]=xlr[i];
x[2*i+1]=xli[i];
}
for(i=0;i<nxl;i++)
{
xm[i]=(int)sqrt((xr[i]*xr[i])+(xi[i]*xi[i]));
}
ChangeOrder(xlr,xli);
radix2(xlr,xli,w1,w2);
for(i=0;i<nxl;i++)
{
int sum=0;
sum=(int)sqrt((xlr[i]*xlr[i])+(xli[i]+xli[i]));
ym[i]=sum;
y[2*i]=(int)sqrt((xlr[i])*(xlr[i]));
y[2*i+1]=(int)sqrt((xli[i]*xli[i]));
}
}
void radix2(int xr[],int xi[],short wr[],short wi[])
{
short L,B,J,P,k;
int rPartKB,iPartKB;
for(L=1;L<=M;L++)
{
B=(int)(pow(2,(L-1))+0.5);
for(J=0;J<=B-1;J++)
{
P=J*((int)(pow(2,(M-L))+0.5));
for(k=J;k<=N-1;k+=(int)(pow(2,L)+0.5))
{
rPartKB=xr[k+B]*wr[P]+xi[k+B]*wi[P];
iPartKB=xi[k+B]*wr[P]-xr[k+B]*wi[P];
rPartKB>>=15;
iPartKB>>=15;
xr[k+B]=xr[k]-rPartKB;
xi[k+B]=xi[k]-iPartKB;
xr[k]=xr[k]+rPartKB;
xi[k]=xi[k]+iPartKB;
}
}
}
}
void ChangeOrder(int xr[],int xi[])
{
int LH,N1,I,J,K;
short 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;
}
}