#include "stdio.h"
#include "math.h"

#define N 21
#define PI 3.1415926
#define fm  35
#define dt 0.004
#define Hf 100

void dft(xr,xi,flag)
float xr[N],xi[N];
int flag;
{ float XR[N],XI[N];
  int k,n;
  float sum1,sum2,cita;
  for(k=0;k<=N-1;k++)
  { sum1=0.0;
    sum2=0.0;
    { for(n=0;n<=N-1;n++)
        {    cita=2.0*PI/N*n*k;
            sum1=sum1+xr[n]*cos(cita)+flag*xi[n]*sin(cita);
            sum2=sum2-flag*xr[n]*sin(cita)+xi[n]*cos(cita);
        }
      XR[k]=sum1;
      XI[k]=sum2;
    }
  }
  if(flag==1)
      for(k=0;k<=N-1;k++)
        {
            xr[k]=XR[k];
            xi[k]=XI[k];
        }
  else
      for(k=0;k<=N-1;k++)
        {
            xr[k]=XR[k]/N;
            xi[k]=XI[k]/N;
        }
}

main()
{ void dft();
  float xr[N],xi[N]={0};
  int i,n,Nf;
  float t,fs,df;
  FILE *fp1,*fp2;
  fs=1.0/dt;
  df=fs/N;
  Nf=Hf/df;

  fp1=fopen("xt1.xls","w");

  fp2=fopen("Xk1.xls","w");

  for(i=-N/2;i<=N/2;i++)
    { t=i*dt;
      xr[i+N/2]=(1.0-2.0*pow(PI*fm*t,2.0))*exp(-pow(PI*fm*t,2.0));
      fprintf(fp1,"%8.3f \t %8.4f\n ",t*1000,xr[i+N/2]);
    }
  fclose(fp1);

  dft(xr,xi,1);


  for(n=0;n<=Nf;n++)
  { xr[n]=dt*sqrt(xr[n]*xr[n]+xi[n]*xi[n]);
    fprintf(fp2,"%8.3f \t %8.3f\n",n*df,xr[n]);
  }
  fclose(fp2);
}