#include "config.h"
#include<math.h>

struct fushu                //fushu結構體表示復數
{ float re;                 //實部
  float im;                 //虛部
};
void fsmul(struct fushu *p,struct fushu a,struct fushu b) //復數乘法
{ p->re=a.re*b.re-a.im*b.im;
  p->im=a.re*b.im+a.im*b.re;
}
void fsadd(struct fushu * p,struct fushu a,struct fushu b) //復數加法
{ p->re=a.re+b.re;
  p->im=a.im+b.im;
}
void fssub(struct fushu * p,struct fushu a,struct fushu b) //復數減法
{ p->re=a.re-b.re;
  p->im=a.im-b.im;
}
/*****************************************/
/*           DIT的基－2FFT算法           */
/* a: 待變換序列指針                     */
/* l: 待變換序列長度                     */
/* n: FFT點數                            */
/*****************************************/
void FFT(struct fushu *a,uint16 l,uint16 n)
{
  uint16 i,j,nv2,nm1,k,le,lei,ip,m,L;
  struct fushu u,w,t;
  float tmp;
  float  pai=atan(1)*4;
  if(l<n)                  //若序列長度不夠FFT變換點數則補零
    for(i=l;i<n;i++)
     {a[i].re=0;a[i].im=0;}
   L=log(n)/log(2);
   nv2=n/2;
   nm1=n-1;
   j=0;
   for(i=0;i<nm1;i++)     //將順序輸入的序列變換成倒位序
    {if(i<j)
      {t.re=a[j].re;t.im=a[j].im;
       a[j].re=a[i].re;a[j].im=a[i].im;
       a[i].re=t.re;a[i].im=t.im;
      }
     k=nv2;
     while(k<=j)
      {j=j-k;k=k/2;}
     j+=k;
    }
   le=1;                  //進行FFT運算
   for(m=1;m<=L;m++)
   {lei=le;
    le=le*2;
    u.re=1;
    u.im=0;
    tmp=pai/lei;
    w.re=cos(tmp);
    w.im=-sin(tmp);
    for(j=0;j<lei;j++)
    {for(i=j;i<n;i+=le)
     {ip=i+lei;
     fsmul(&t,a[ip],u);
     fssub(&(a[ip]),a[i],t);
     fsadd(&(a[i]),a[i],t);
     }
     fsmul(&u,u,w);
    }
   }
}