#include "math.h"
void fft(double x[],double y[],int n,int sign)
//int n,sign;
//double x[],y[];
{
int i,j,k,l,m,n1,n2;
double c,c1,e,s,s1,t,tr,ti;
for(j=1,i=1;i<16;i++)
{
m=i;
j=2*j;
if(j==n) break;
}
n1=n-1;
for(j=0,i=0;i<n1;i++)
{
if(i<j)
{
tr=x[j];
ti=y[j];
x[j]=x[i];
y[j]=y[i];
x[i]=tr;
y[i]=ti;
}
k=n/2;
while(k<(j+1))
{
j=j-k;
k=k/2;
}
j=j+k;
}
n1=1;
for(l=1;l<=m;l++)
{
n1=2*n1;
n2=n1/2;
e=3.14159265359/n2;
c=1.0;
s=0.0;
c1=cos(e);
s1=-sign*sin(e);
for(j=0;j<n2;j++)
{
for(i=j;i<n;i+=n1)
{
k=i+n2;
tr=c*x[k]-s*y[k];
ti=c*y[k]+s*x[k];
x[k]=x[i]-tr;
y[k]=y[i]-ti;
x[i]=x[i]+tr;
y[i]=y[i]+ti;
}
t=c;
c=c*c1-s*s1;
s=t*s1+s*c1;
}
}
if(sign==-1)
{
for(i=0;i<n;i++)
{
x[i]/=n;
y[i]/=n;
}
}
}
#include "math.h"
#include "stdio.h"
//#include "fft.c"
void main()
{
int i,j,n;
double a1,a2,c,c1,c2,d1,d2,q1,q2,w,w1,w2;
double x[32],y[32],a[32],b[32];
n=32;
a1=0.9;
a2=0.3;
x[0]=1.0;
y[0]=0.0;
for(i=1;i<n;i++)
{
x[i]=a1*x[i-1]-a2*y[i-1];
y[i]=a2*x[i-1]+a1*y[i-1];
}
printf("\nCOMPLEX INPUT SEQUENCE\n");
for(i=0;i<n/2;i++)
{
for(j=0;j<2;j++)
printf("%10.7f+J%10.7f",x[2*i+j],y[2*i+j]);
printf("\n");
}
q1=x[n-1];
q2=y[n-1];
for(i=0;i<n;i++)
{
w=6.28318530718/n*i;
w1=cos(w);
w2=-sin(w);
c1=1.0-a1*w1+a2*w2;
c2=a1*w2+a2*w1;
c=c1*c1+c2*c2;
d1=1.0-a1*q1+a2*q2;
d2=a1*q2+a2*q1;
a[i]=(c1*d1+c2*d2)/c;
b[i]=(c2*d1-c1*d2)/c;
}
printf("\n THEORETICAL DFT\n");
for(i=0;i<n/2;i++)
{
for(j=0;j<2;j++)
printf("%10.7f+J%10.7f",a[2*i+j],b[2*i+j]);
printf("\n");
}
fft(x,y,n,1);
printf("\nDISCRETE FORIER TRANSFORM\n");
for(i=0;i<n/2;i++)
{
for(j=0;j<2;j++)
{
printf("%10.7f+J%10.7f",x[2*i+j],y[2*i+j]);
}
printf("\n");
}
fft(x,y,n,-1);
printf("\n INVERSE DISCRETE FOURIER TRANSFORM\n");
for(i=0;i<n/2;i++)
{
for(j=0;j<2;j++)
printf("%10.7f+J%10.7f",x[2*i+j],y[2*i+j]);
printf("\n");
}
}
