#include <stdio.h>
#include <math.h>
#define pi 3.1415926
#define SIG_LEN 256
typedef struct{
float real;
float image;
}Complex;
Complex ComMul(Complex p1,Complex p2);
Complex ComAdd(Complex p1,Complex p2);
Complex ComSub(Complex p1,Complex p2);
int Reform(float **in,int len);
int BitRevers(int src,int size);
void Displace(float *in,int size,int M);
void CalcW(Complex *w,int len);
Complex *fft_1d(float *in,int len);
main()
{
int slen,i;
FILE *fp,*fout;
float *data;
Complex *out;
data = malloc(SIG_LEN*sizeof(float));
if((fp=fopen("e:\input","r"))==NULL)
{
printf("cannot open file");
exit(0);
}
if((fout=fopen("e:\output","w+"))==NULL)
{
printf("cannot open file");
exit(0);
}
for(i=0;i<SIG_LEN;i++)
fscanf(fp,"%f\n",&data[i]);
out = fft_1d(data,SIG_LEN);
for(i=0;i<SIG_LEN;i++)
fprintf(fout,"%f\n",sqrt(pow(out[i].real,2)+pow(out[i].image,2)));
}
Complex ComMul(Complex p1,Complex p2)
{
Complex res;
res.real = p1.real*p2.real - p1.image*p2.image;
res.image = p1.real*p2.image + p1.image*p2.real;
return res;
}
Complex ComAdd(Complex p1,Complex p2)
{
Complex res;
res.real = p1.real + p2.real;
res.image = p1.image + p2.image;
return res;
}
Complex ComSub(Complex p1,Complex p2)
{
Complex res;
res.real = p1.real - p2.real;
res.image = p1.image - p2.image;
return res;
}
int BitRevers(int src,int size)
{
int temp=src;
int dst=0;
int i=0;
for(i=size-1;i>=0;i--)
{
dst=((temp&0x1)<<i)|dst;
temp=temp>>1;
}
return dst;
}
void Displace(float *in,int size,int M)
{
int i;
int new_i;
float t;
FILE *fp;
for(i=1;i<size;i++)
{
new_i=BitRevers(i,M);
if(new_i>i)
{
t=in[i];
in[i]=in[new_i];
in[new_i]=t;
}
}
if((fp=fopen("e:\displace","w+"))==NULL)
{
printf("cannot open file");
exit(0);
}
for(i=0;i<size;i++)
fprintf(fp,"%f\n",in[i]);
}
void CalcW(Complex *w,int N)
{
int i;
FILE *fp;
for(i=0;i<N;i++)
{
w[i].real = cos(pi*i/N);
w[i].image = -1*sin(pi*i/N);
}
if((fp=fopen("e:\w","w+"))==NULL)
{
printf("cannot open file");
exit(0);
}
for(i=0;i<N;i++)
{fprintf(fp," .word %ld\n",(long)(w[i].real*16384));
fprintf(fp," .word %ld\n",(long)(w[i].image*16384));
}
}
int Reform(float **in,int len)
{
int i=0;
int w=1;
while(w*2<=len) {w=w*2;i++;}
if(w<len)
{
*in = realloc(*in,w*2*sizeof(float));
for(i=len;i<w*2;i++)
in[i]=0;
return i+1;
}
return i;
}
Complex *fft_1d(float *data,int len)
{
int M,i=0,j,k;
int GroupNum;
int CellNum;
int reallen;
int pos1,pos2;
float *in;
Complex *res,*t;
Complex *w,mul;
in = data;
M = Reform(&in,len); /*let the data length is pow of 2*/
reallen = pow(2,M);
res = malloc(reallen*sizeof(Complex));
w = malloc(reallen*sizeof(Complex)/2);
CalcW(w,reallen/2);
Displace(in,reallen,M);
while(i<reallen) {res[i].real = in[i];res[i].image=0;i++;}
GroupNum = reallen/2;
CellNum = 1;
for(i=0;i<M;i++)
{
for(j=0;j<GroupNum;j++)
{
for(k=0;k<CellNum;k++)
{
pos1 = j*CellNum*2 + k;
pos2 = pos1 + CellNum;
mul = ComMul(res[pos2],w[k*GroupNum]);
res[pos2] = ComSub(res[pos1] ,mul);
res[pos1] = ComAdd(res[pos1] ,mul);
}
}
GroupNum = GroupNum/2;
CellNum = CellNum*2;
}
return res;
}