//---------------------------------------------------------------------------
/*
pr-----实型一维数组，长度为n。当n=0时，存放n个采样输入的实部，返回时存放离散傅立叶变换的模；当l＝1 时，存放傅立叶变换的n个实部，返回时存放逆傅立叶变换的模。 pi-----实型一维数组，长度为n。当n=0时，存放n个采样输入的虚部，返回时存放离散傅立叶变换的幅角；当l＝ 1时，存放傅立叶变换的n个虚部，返回时存放逆傅立叶变换的幅角。其中幅角的单位为度。 n------整形。输入的点数。 k------整形。满足n=2~k（2的k次方）。 fr-----实型一维数组，长度为n。当l=0时，返回傅立叶变换的实部；当l＝1时，返回逆傅立叶变换的实部。 fi-----实型一维数组，长度为n。当l=0时，返回傅立叶变换的虚部；当l＝1时，返回逆傅立叶变换的虚部。 l------整形。若l＝0，表示要求本函数计算傅立叶变换；若l＝1，表示要求本函数计算逆傅立叶变换。
*/
void kbfft(double *pr,double *pi,int n,int k,double *fr,double *fi,int l)
{
 int it,m,is,i,j,nv,l0;
 double p,q,s,vr,vi,poddr,poddi;
 
 for(it=0;it<=n-1;it++)
 {
  m=it;
  is=0;
  for(i=0;i<=k-1;i++)
  {
   j=m/2;
   is=2*is+(m-2*j);
   m=j;
  }
  if(l==1)
                {
                fr[it]=pr[is]*n/4;
  fi[it]=pi[is]*n/4;
                }
                else
                {
                fr[it]=pr[is]/n*4;
  fi[it]=pi[is]/n*4;
                }  }  pr[0]=1.0;
 pi[0]=0.0;
 p=6.283185306/(1.0*n);
 pr[1]=cos(p);
 pi[1]=-sin(p);
 
 if(l!=0)
  pi[1]=-pi[1];
 
 for(i=2;i {
  p=pr[i-1]*pr[1];
  q=pi[i-1]*pi[1];
  s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
  pr[i]=p-q;pi[i]=s-p-q;
 }
 
 for(it=0;it<=n-2;it=it+2)
 {
  vr=fr[it];
  vi=fi[it];
  fr[it]=vr+fr[it+1];
  fi[it]=vi+fi[it+1];
  fr[it+1]=vr-fr[it+1];
  fi[it+1]=vi-fi[it+1];
 }
 
 m=n/2;
 nv=2;
 
 for(l0=k-2;l0>=0;l0--)
 {
  m=m/2;
  nv=2*nv;
  for(it=0;it<=(m-1)*nv;it=it+nv)
   for(j=0;j<=(nv/2)-1;j++)
   {
    p=pr[m*j]*fr[it+j+nv/2];
    q=pi[m*j]*fi[it+j+nv/2];
    s=pr[m*j]+pi[m*j];
    s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
    poddr=p-q;
    poddi=s-p-q;
    fr[it+j+nv/2]=fr[it+j]-poddr;
    fi[it+j+nv/2]=fi[it+j]-poddi;
    fr[it+j]=fr[it+j]+poddr;
    fi[it+j]=fi[it+j]+poddi;
   }
 }
 
 if(l!=0)
 {
  for(i=0;i<=n-1;i++)
  {
   fr[i]=fr[i]/(1.0*n);
   fi[i]=fi[i]/(1.0*n);
  }
 }
 return;
}
//--------------------