/*****************fft programe*********************/
/*混合基FFT算法*/
//by SC07023029 呂嫻娜 
#include "math.h"
#include "stdio.h"

struct compx 
{ double real;
  double imag;
} compx ;


struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}
struct compx add(struct compx a,struct compx b)
{struct compx c;
 c.real=a.real+b.real;
 c.imag=a.imag+b.imag;
return(c);
}

void FFT(struct compx *xin,int N,int r1,int r2 )

{
int k0,k1,n0,n1;
double p ,ps ;
float pi;
struct compx v[N],w,t;

{

for (k1=0;k1<=r1-1;k1++)
    {for(k0=0;k0<=r2-1;k0++)
        {  
         pi=3.14159;
         for(n0=0;n0<=r2-1;n0++)
         {
             for(n1=0;n1<=r1-1;n1++)
             {p=r2*n1*k0;
              ps=2*pi/N*p;
              w.real=cos(ps);
              w.imag=-sin(ps);
              t=EE(xin[n1*r2+n0],w);
              v[n0]=add(v[n0],t);
              }
              p=n0*k0;
              ps=2*pi/N*p;
              w.real=cos(ps);
              w.imag=-sin(ps);
              v[n0]=EE(v[n0],w);
              p=r1*n0*k1;
              ps=2*pi/N*p;
              w.real=cos(ps);
              w.imag=-sin(ps);
         v[n0]=EE(v[n0],w);   
         xin[k1*r1+k0]=add(xin[k1*r1+k0],v[n0]) ;
         }
       }
    }
  }  
return ;
} 


/*****************main programe********************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

float  result[257];
struct  compx s[257];
int   Num=16;
int r1=4,r2=4;
const float pp=3.14159;

main()
{

int i;
printf("原序列：\n");
for(i=0;i<16;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0;
printf("%.4f    ",s[i].real);
}

FFT(s,Num,r1,r2);

for(i=0;i<16;i++)
{
printf("%.4f",s[i].real);
printf("+%.4fj    ",s[i].imag);
}
printf("\n功率密度:\n");
for(i=0;i<16;i++)
{
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
printf("+%.4f    ",result[i]);
}
scanf("2");

}