float FFT_power_spectrum(void)
{ 
  /*
  FFT(dat);                                      //FFT計算
  //////////////////////////////////////////////////////////////
  for(i=0;i<(NNN/2);i++)
  Power_spectrum_dat[i+1]=sqrt(dat[i*2+1]*dat[i*2+1]+dat[i*2+2]*dat[i*2+2]);
  Power_spectrum_dat[0]=Power_spectrum_dat[2];
  K=2;
  for(i=3;i<=(NNN/2);i++) 
    {
       if(Power_spectrum_dat[i]>Power_spectrum_dat[0]) 
          {
             Power_spectrum_dat[0]=Power_spectrum_dat[i];
             K=i;
          }
       else continue;
    }
////////////////////////////////////////////////
  i=K; 

  theta1=0,theta2=0,theta3=0;                    //基于相角判據的修正Rife算法
  theta2=atan(dat[i*2]/dat[i*2-1]);
  theta3=atan(dat[i*2+2]/dat[i*2+1]);
  theta1=atan(dat[i*2-2]/dat[i*2-3]);
  if(fabsf(theta2-theta1)>fabsf(theta2-theta3))      
  {
    K+=(Power_spectrum_dat[i+1]/(Power_spectrum_dat[i]+Power_spectrum_dat[i+1]));
  }
  else if(fabsf(theta2-theta1)==fabsf(theta2-theta3))       
  {
    K+=0;
  }
  else
  {
    K-=(Power_spectrum_dat[i-1]/(Power_spectrum_dat[i]+Power_spectrum_dat[i-1]));
  }
  */
/*  
  K+=(-Power_spectrum_dat[i-1]+Power_spectrum_dat[i+1])/(2*Power_spectrum_dat[i]-Power_spectrum_dat[i-1]+Power_spectrum_dat[i+1]);    //拋物線內插法
*/ 
/*  
  if(Power_spectrum_dat[i-1]<=Power_spectrum_dat[i+1])     //雙線性對稱頻率內插??
  {
    K+=(Power_spectrum_dat[i+1]-Power_spectrum_dat[i-1])/(Power_spectrum_dat[i]-Power_spectrum_dat[i-1]);
  }
  else
  {
    K-=(Power_spectrum_dat[i+1]+Power_spectrum_dat[i-1])/(Power_spectrum_dat[i]-Power_spectrum_dat[i+1]);
  }
*/  
/*
  float delta1,delta2;                           //DFT變換系數的實部構造頻率修正項的算法??
  delta1=0,delta2=0;
  delta1=dat[i*2+1]/(dat[i*2+1]-dat[i*2-1]);
  delta2=dat[i*2-3]/(dat[i*2-1]-dat[i*2-3]);
  if(delta1>0&delta2>0)
  {
    K+=delta1;
  }
  else if(dat[i*2+1]==dat[i*2-3])
  {
    K+=0;
  }
  else
  {
    K+=delta2;
  }
*/
/*  
  K=(Power_spectrum_dat[i-2]*(i-2)+Power_spectrum_dat[i-1]*(i-1)+Power_spectrum_dat[i-0]*(i-0)+Power_spectrum_dat[i+1]*(i+1)+Power_spectrum_dat[i+2]*(i+2))/
  (Power_spectrum_dat[i-2]+Power_spectrum_dat[i-1]+Power_spectrum_dat[i]+Power_spectrum_dat[i+1]+Power_spectrum_dat[i+2]);        //能量重心法5點
*/
/* 
  float theta1,theta2,deltak;
  theta1=0,theta2=0,deltak=0;
  for(i=3;i<=(NNN/4);i++)                        //相位差法頻譜校正1??
  {
    if(Power_spectrum_dat[i]>Power_spectrum_dat[0]) 
    {
      Power_spectrum_dat[0]=Power_spectrum_dat[i];
    }
    else continue;
  }
  theta1=atan2(dat[i*2],dat[i*2-1]);
  for(j=NNN/4+1;j<=(NNN/2);j++) 
  {
    if(Power_spectrum_dat[j]>Power_spectrum_dat[NNN/4+1]) 
    {
      Power_spectrum_dat[NNN/4+1]=Power_spectrum_dat[j];
    }
    else continue;
  }
  theta2=2(dat[j*2],dat[j*2-1]);
  deltak=(theta1-theta2)/(2*M_PI);
  if(deltak<-0.5)
  {
    deltak=deltak+1;
  }
  else if(deltak>0.5)
  {
    deltak=deltak-1;
  }
  else
  {
    deltak=deltak;    
  }
  K+=deltak;
*/
/*
  unsigned int km0,km1;
  km0=0,km1=0;
  for(i=3;i<=(NNN/2);i++)                        //相位差法頻譜校正2??
  {
    if(Power_spectrum_dat[i]>Power_spectrum_dat[0]) 
    {
      Power_spectrum_dat[0]=Power_spectrum_dat[i];
      km0=i;
    }
    else continue;
  }
  for(j=3;j<=(NNN/0);j++) 
  {
    if(Power_spectrum_dat[j]>Power_spectrum_dat[NNN/4+1]) 
    {
      Power_spectrum_dat[NNN/4+1]=Power_spectrum_dat[j];
      km1=j;
    }
    else continue;
  }
  f_fft=((fs*fs_edit)/(NNN/4))*(((NNN/2-1)/(NNN/2))*km0-((NNN/4-1)/(NNN/4))*km1);
*/
/*
  if(Power_spectrum_dat[i-1]<=Power_spectrum_dat[i+1])     //Rife
  {
    K+=Power_spectrum_dat[i+1]/(Power_spectrum_dat[i]+Power_spectrum_dat[i+1]);
  }
  else
  {
    K-=Power_spectrum_dat[i-1]/(Power_spectrum_dat[i]+Power_spectrum_dat[i-1]);
  }
*/
/////////////////////////////////////////////////
 
  /*f_fft=(float)(((K-1)*fs)/NNN);  
  return f_fft;
}
void FFT(float xxx[])
{   
  n=NNN; q=1;
  for(p=1;p<n;p+=2)
  {
    if(q>p)
    {
      SWAP(xxx[q],xxx[p]);
      SWAP(xxx[q+1],xxx[p+1]);
    }
    m=n>>1;
    while(m>=2&&q>m)
    {
      q-=m;
      m>>=1;
    }
    q+=m;
  }
  mmax=2;
  while(n>mmax)
  {
    istep=2*mmax;
    theta=2.0*M_PI/mmax;
//////////////////////////////
    wtemp=sin((float)(0.5*theta));
    wpr=-2.0*wtemp*wtemp;  
//////////////////////////////          
    wpi=sin((float)theta);
    for(wr=1.0,wi=0.0,m=1;m<mmax;m+=2)
    {
      for(p=m;p<=n;p+=istep)
      {
        q=p+mmax;
        tempr=wr*xxx[q]-wi*xxx[q+1];
        tempi=wr*xxx[q+1]+wi*xxx[q];
        xxx[q]=xxx[p]-tempr;
        xxx[q+1]=xxx[p+1]-tempi;
        xxx[p]+=tempr;
        xxx[p+1]+=tempi;
      }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
    }
    mmax=istep;
  }
  */
}