#include         <stdio.h>
%A #include         <math.h>
%A
%A #define PI 3.14159265358979
%A #define numlength 1024*64
%A #define TIME          1000
%A
%A void        radix2(int n, short *xy, short *w);
%A void        bitrev_cplx(int *x, short *index, int nx);
%A void        digitrev_index(short *index, int n, int radix);
%A void        dif_2(void);
%A void FFT(float xr[],float xi[],int MA2);
%A
%A
%A short index1[64],w1[numlength];
%A short x1[numlength*2];
%A int x2[numlength];
%A int radix;
%A int nx1;
%A int i;
%A float rm[1024],im[1024];
%A float rm1[512],im1[512];
%A         
%A /*======================主程序====================================*/
%A
%A void         main()
%A {        
%A         int cnt;        
%A         *(int *)0x10070 = 0x12345678;
%A         *(int *)0x10100 = 0;                
%A         dif_2();
%A         for(i=0;i<512;i++)
%A         {
%A                 rm[i] = cos(PI*i*5)*128.0;
%A                 im[i] = 0;
%A         }
%A         
%A         for (cnt=0;cnt<1000;cnt++)
%A         {
%A                 for(i=0;i<512;i++)
%A                 {
%A                         rm1[i] = rm[i];
%A                         im1[i] = im[i];
%A                 }        
%A                 FFT(rm1,im1,512);
%A         }
%A         *(int *)0x10100 = 0x55555555;
%A
%A }
%A
%A void FFT(float xr[],float xi[],int MA2)
%A {
%A         int invr=1;
%A         float ur,ui,tr,ti;
%A         int i,j,k,l,m,le,lei;
%A
%A         j=0;
%A         for(i=0;i<MA2-1;i++)
%A         {
%A                 if(i<j)
%A                 {
%A                 tr=xr[i];ti=xi[i];
%A                 xr[i]=xr[j];xi[i]=xi[j];
%A                 xr[j]=tr;xi[j]=ti;
%A                 }
%A                 k=MA2>>1;
%A                 while(j>=k)
%A                 {
%A                         j-=k;
%A                         k>>=1;
%A                 }
%A            j+=k;
%A         }
%A
%A         i=MA2;
%A         for(l=0;i>1;l++) i>>=1;
%A         le=1;
%A         for(i=0;i<l;i++)
%A         {
%A                 lei=le;///2;
%A                 le<<=1;//2*le;
%A                 for(j=0;j<lei;j++)
%A                 {
%A                         float dtemp = j*PI/lei;
%A                         ur=cos(dtemp);//invr*
%A                         ui=-sin(dtemp);//invr*
%A                         for(k=j;k<MA2;k+=le)
%A                         {
%A                                 m=k+lei;
%A                                 tr=xr[m]*ur-xi[m]*ui;
%A                                 ti=xr[m]*ui+xi[m]*ur;
%A                                 xr[m]=xr[k]-tr;
%A                                 xi[m]=xi[k]-ti;
%A                                 xr[k]+=tr;
%A                                 xi[k]+=ti;
%A                         }
%A                 }
%A         }
%A         xi[0]=0;
%A         xr[0]=0;
%A //        xi[1]=0;
%A //        xr[1]=0;
%A }
%A
%A void        dif_2(void)
%A {
%A         float delta;
%A         int cnt;
%A         
%A         *(int *)0x10070 = 0x12345678;
%A         *(int *)0x10100 = 0;
%A         
%A         radix = 2;                        // 基数 --- 基2FFT
%A         nx1 = numlength;                //FFT计算点数
%A
%A         //构造要做FFT的序列，按实部0，虚部0，实部1，虚部1。。。顺序构造。
%A         //本例中实部为三个频率的信号。虚部为0
%A         
%A         for(i=0;i<nx1;i++)
%A         {
%A                 x1[2*i] = (short)((cos(PI*i/10.0) + cos(PI*i/20.0) + cos(PI*i/50.0))*0x80);
%A                 x1[2*i+1] = 0;
%A         }
%A         
%A         delta = 2*PI/nx1;
%A         for(i=0;i<nx1/2;i++)
%A         {
%A                 w1[2*i] = 32767*(-cos((float)i*delta));
%A                 w1[2*i+1] = 32767*(-sin((float)i*delta));
%A         }        
%A
%A         for(cnt=0;cnt<1000;cnt++)        
%A         {
%A                 //init the index for FFT bitrev
%A                 digitrev_index(index1,nx1,radix);
%A         
%A                 radix2(nx1, x1, w1);
%A                         
%A                 for(i=0;i<nx1;i++)
%A                 {
%A                         x2[i] = sqrt(x1[2*i]*x1[2*i] + x1[2*i+1]*x1[2*i+1]);
%A                 }
%A                         
%A                 //倒序
%A                         
%A                 bitrev_cplx(x2, index1, nx1);
%A         }
%A
%A         *(int *)0x10100 = 0x55555555;
%A }  
%A
%A //生成倒序表子程序
%A
%A void        digitrev_index(short *index, int n, int radix)
%A {
%A         int i,j,k;
%A         short nbits, nbot, ntop, ndiff, n2, raddiv2;
%A         
%A         nbits = 0;
%A         i = n;
%A         while(i > 1)
%A         {
%A         i = i >> 1;
%A         nbits++;
%A         }
%A         raddiv2 = radix >> 1;
%A         nbot = nbits >> raddiv2;
%A         nbot  = nbot << raddiv2 - 1;
%A         ndiff = nbits&raddiv2;
%A         ntop = nbot + ndiff;
%A         n2 = 1<<ntop;
%A         index[0] = 0;
%A         for(i=1,j=n2/radix+1;i<n2-1;i++)
%A         {
%A                 index[i] = j - 1;
%A                 for(k=n2/radix;k*(radix-1)<j;k/=radix)
%A                 {
%A                 j -= k*(radix - 1);
%A                 }
%A                 j += k;
%A         }
%A         index[n2-1] = n2 - 1;
%A }
%A
%A //倒序子程序
%A void        bitrev_cplx(int *x, short *index, int nx)
%A {
%A         int i;
%A         short i0, i1, i3;
%A         short j0, j1, j3;
%A         int xi0, xi1, xi3;
%A         int xj0, xj1, xj3;
%A         short t;
%A         int a, b, ia, ib, ibs;
%A         int mask;
%A         int nbits, nbot, ntop, ndiff, n2, halfn;
%A         
%A         nbits = 0;
%A         i = nx;
%A         
%A         while(i > 1)
%A         {
%A                 i = i >> 1;
%A                 nbits++;
%A         }
%A         
%A         nbot = nbits >> 1;
%A         ndiff = nbits&1;
%A         ntop = nbot + ndiff;
%A         n2 = 1 << ntop;
%A         mask = n2 - 1;
%A         halfn = nx >> 1;
%A         
%A         for(i0=0;i0<halfn;i0+=2)
%A         {
%A                 b = i0&mask;
%A                 a = i0 >> nbot;
%A                 
%A                 if(!b)
%A                 {
%A                         ia = index[a];
%A                 }
%A                 ib = index[b];
%A                 ibs = ib << nbot;
%A                 
%A                 j0 = ibs + ia;
%A                 
%A                 t = i0 < j0;
%A                 xi0 = x[i0];
%A                 xj0 = x[j0];
%A                 
%A                 if(t)
%A                 {
%A                         x[i0] = xj0;
%A                         x[j0] = xi0;
%A                 }
%A                 
%A                 i1 = i0 + 1;
%A                 j1 = j0 + halfn;
%A                 xi1 = x[i1];
%A                 xj1 = x[j1];
%A                 x[i1] = xj1;
%A                 x[j1] = xi1;
%A                 
%A                 i3 = i1 + halfn;
%A                 j3 = j1 + 1;
%A                 xi3 = x[i3];
%A                 xj3 = x[j3];
%A                 
%A                 if(t)
%A                 {
%A                         x[i3] = xj3;
%A                         x[j3] = xi3;
%A                 }
%A         }
%A }
%A
%A //频率抽取FFT算法函数
%A //xy[] ------ input and output sequences(dim-n) (input/output)
%A // n   ------ FFT size (input)
%A // w[] ------ FFT coefficients(dim-n/2) (input)
%A
%A //DESCRIPTION
%A // This routine is used to compute FFT of a complex sequence
%A // of size n, a power of 2, with "decimation-in-frequency
%A // decomposition"method,ie,the output is in bit-reversed
%A // order.Each complex value is with interleaved 16-bit real
%A // and imaginary parts.
%A
%A void        radix2(int n, short xy[], short w[])
%A {
%A         short n1,n2,ie,ia,i,j,k,l;
%A         short xt,yt,c,s;
%A         
%A         n2 = n;
%A         ie = 1;
%A         
%A         for(k=n;k>1;k=(k>>1))
%A         {
%A                 n1 = n2;
%A                 n2 = n2 >> 1;
%A                 ia = 0;
%A                 
%A                 for(j=0;j<n2;j++)
%A                 {
%A                         c = w[2*ia];
%A                         s = w[2*ia+1];
%A                         ia = ia + ie;
%A                         
%A                         for(i=j;i<n;i+=n1)
%A                         {
%A                                 l = i + n2;
%A                                 xt = xy[2*l] - xy[2*i];
%A                                 xy[2*i] = xy[2*i] + xy[2*l];
%A                                 yt = xy[2*l+1] -xy[2*i+1];
%A                                 xy[2*i+1] = xy[2*i+1] + xy[2*l+1];
%A                                 xy[2*l] = (c*xt + s*yt) >> 15;
%A                                 xy[2*l+1] = (c*yt - s*xt) >> 15;
%A                         }
%A                 }
%A                 ie = ie << 1;
%A         }
%A }
%A
%A
%A
%A
%A
%A
%A
%A
%A
%A
%A
%A
%A
%A
%A %F