/**********************************************************************
/*文件名：fft.cpp
/* 
/*   fft源文件
/*   修改時間:2004.9.2
/*   作者：* * *
/**********************************************************************/
#include "fft.h"

void splfft( double xr[],double xi[], int n,int isign)
{/*----------------------------------------------------------------------
 功能：實現FFT或其逆變換

 輸入參數說明：
 xr[] :待處理數據的實部數組,xr[0]----xr[n-1]
 xi[] :待處理數據的虛部數組,xi[0]----xi[n-1]
 n    :數組大小
 isign:如果 isign=-1 實現FFT正變換，并把變換結果的實部保存到xr[]，虛部保存到xi[]
       如果 isign=+1 實現FFT逆變換，并把變換結果的實部保存到xr[]，虛部保存到xi[]
 
 注意：n必須是2的正整數次冪，如：256、512、1024
 -------------------------------------------------------------------------*/
    double xtr,xti;
        double es,e,a,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3;
        int m,n2,k,n4,j,is,id,i1,i2,i3,i0,n1,i,nn;
        for(m=1;m<=16;m++)
    	  {
    	    nn=(int)pow(2,m);
    	    if(n==nn)break;
    	  }
	if(m>16)
	  {
        printf(" N is not a power of 2 ! \n");
	    return;
	  }
	    n2=n*2;
        es=-isign*atan(1.0)*8.0;
   for(k=1;k<m;k++)
     {
            n2=n2/2;
            n4=n2/4;
            e=es/n2;
            a=0.0;
            for(j=0;j<n4;j++)
               {
                a3=3*a;
                cc1=cos(a);
                ss1=sin(a);
                cc3=cos(a3);
                ss3=sin(a3);
                a=(j+1)*e;
                is=j;
                id=2*n2;
                 do
                  {
                    for(i0=is;i0<n;i0+=id)
                       {
                           i1=i0+n4;
                            i2=i1+n4;
                            i3=i2+n4;
                            r1=xr[i0]-xr[i2];
                            s1=xi[i0]-xi[i2];
                            r2=xr[i1]-xr[i3];
                            s2=xi[i1]-xi[i3];
                            xr[i0]+=xr[i2];xi[i0]+=xi[i2];
                            xr[i1]+=xr[i3];xi[i1]+=xi[i3];
                		    if(isign!=1)
                    		    {
                    		      s3=r1-s2;
                    		      r1+=s2;
                    		      s2=r2-s1;
                    		      r2+=s1;
                    		    }
                    		 else
                    		    {
                    		      s3=r1+s2;
                    		      r1=r1-s2;
                    		      s2=-r2-s1;
                    		      r2=-r2+s1;
                                 }
                    		xr[i2]=r1*cc1-s2*ss1;
                    		xi[i2]=-s2*cc1-r1*ss1;
                    		xr[i3]=s3*cc3+r2*ss3;
                    		xi[i3]=r2*cc3-s3*ss3;
                        }
                   is=2*id-n2+j;
                  id=4*id;
		  }while(is<n-1);
      }
  }
/*   ------------ special last stage -------------------------*/
        is=0;
        id=4;
      do
      {
        for(i0=is;i0<n;i0+=id)
           {i1=i0+1;
            xtr=xr[i0];
            xti=xi[i0];
            xr[i0]=xtr+xr[i1];
            xi[i0]=xti+xi[i1];
            xr[i1]=xtr-xr[i1];
            xi[i1]=xti-xi[i1];
            }
        is=2*id-2;
        id=4*id;
       }while(is<n-1);
        j=1;
        n1=n-1;
        for(i=1;i<=n1;i++)
	   {
    	   if(i<j)
    	       {
        		 xtr=xr[j-1];
        		 xti=xi[j-1];
        		 xr[j-1]=xr[i-1];
        		 xi[j-1]=xi[i-1];
        		 xr[i-1]=xtr;
        		 xi[i-1]=xti;
        		 }
    	      k=n/2;
            do
        	  {
            		if(k>=j)
                		 break;
            	    j-=k;
            	    k/=2;
        	    }while(1);
            j+=k;
        }
        if(isign==-1)
           return;
        for(i=0;i<n;i++)
    	   {
    	    xr[i]/=(double)n;
            xi[i]/=(double)(n);
            }
 }