// fft5.cpp : Defines the entry point for the console application.
//
//本程序是為了驗證諧波的算法.
#include "stdafx.h"//#include <complex.h>
#include <math.h>
#include <iostream.h>
//#define N 256
void FFT(struct compx *,int);
struct compx
{
 float real,imag;
};
struct compx s[257];
struct compx xout[257];
struct compx EE(struct compx,struct 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);
}
float pi;
void main()
{
 int times;
 float pi=3.1415926535;
 
 float y;
 
 float y1;
 float num;
 float num1;
 float Vma;
 float Vmb;
 float Vm;
 float z;
 
 const float timeT=2*pi/256;
 
// num=0.0;
//    num1=0.0;
 for (times=0;times<257;times++)
 {
  //s[times].real=times;
  //s[times].imag=times;
  //s[times].real=cos(2*3.1415926*times/256);
  //s[times].imag=-sin(2*3.1415926*times/256);
 }
 for (times=1;times<257;times++)
 {
    s[times].real=(float)(sin((times-1)*timeT)+0.3*sin(5*(times-1)*timeT)+0.01*sin(9*(times-1)*timeT))*sin(pi/3);//*sin(5*times*timeT);
    s[times].imag=0.0;//(float)(sin((times-1)*timeT)+0.3*sin(5*(times-1)*timeT)+0.01*sin(9*(times-1)*timeT))*cos(pi/3);//*cos(5*times*timeT);
   
    //y=(sin((2*pi*times)/256)+0.3*sin((2*times*pi*5)/(256)))*sin((2*times*pi*5)/(256));
    //num=y+num;
    ///num1=num1+y1;
 }
 //Vma=2*num/N;
 //Vmb=2*num1/N;
 //Vm=Vma*Vma+Vmb*Vmb;
    //Vm=sqrt(Vm);
 //z=Vm;
 //
 //s[256].real=0.0;
 //s[256].imag=0.0;
 FFT(s,256);
    for (times=0;times<257;times++)
 {
  y=s[times].real*2/256;
  y1=s[times].imag*2/256;
  y=y*y;
  y1=y1*y1;
  num=y+y1;
  xout[times].real=sqrt(num);
  xout[times].imag=0.0;
 }
 while (1);
}

void FFT(struct compx *xin,int N)
{
 long HighStart,LowStart,HighEnd,LowEnd;
    long numhigh,numlow;

 //const double pi=3.141592653589793;
    int f,m,nv2,nm1,i,k,j=1,l;
 struct compx v,w,t;
 __asm             
 {
  RDTSC
  mov HighStart, edx
  mov LowStart, eax
 }
 nv2=N/2;
 f=N;
 for (m=1;(f=f/2)!=1;m++)
 {
  ;
 }
 nm1=N-1;
 for (i=1;i<=nm1;i++)
 {
  if (i<j)
        {
   t.imag=xin[j].imag;
   t.real=xin[j].real;
   xin[j].imag=xin[i].imag;
   xin[j].real=xin[i].real;
   xin[i].imag=t.imag;
   xin[i].real=t.real;
  }
  k=nv2;
  while (k<j)
  {
   j=j-k;k=k/2;
  }
  j=j+k;
 }
 for (i=0;i<257;i++)
 {
        xout[i]=xin[i];
 }
 //fft
 {
  int le,lei,ip;
  float pi;//,x,y;
  for (l=1;l<=m;l++)
  {
   le=int(pow(2,l));
   lei=le/2;
   pi=3.141592653;
      v.real=1.0;
   v.imag=0.0;
   w.real=(float)cos(pi/lei);
   w.imag=(float)-sin(pi/lei);
   for (j=1;j<=lei;j++)
   {
    for (i=j;i<=N;i=i+le)
    {
     ip=i+lei;
     t=EE(xin[ip],v);
     xin[ip].real=xin[i].real-t.real;
     xin[ip].imag=xin[i].imag-t.imag;
     xin[i].real=xin[i].real+t.real;
     xin[i].imag=xin[i].imag+t.imag;
    }
    v=EE(v,w);
   }
  }
 }
__asm
 {
  RDTSC
  mov HighEnd, edx
  Mov LowEnd,  eax
  ;獲取兩次計數器值得差
  sub eax,  LowStart
  cmp    eax,  0       ; 如果低32的差為負則求返,因為第二次取得永遠比第一次的大
  jg     L1
  neg     eax
  jmp     L2
            L1: mov numlow,  eax
            L2: sbb edx,  HighStart
  mov numhigh, edx
 
 }
        //把兩個計數器值之差放在一個64位的整形變量中
        //先把高32位左移32位放在64的整形變量中,然后再加上低32位
 __int64  timer =(numhigh<<32) + numlow;
         //輸出代碼段運行的時鐘周期數
         //以頻率1.1Gcpu為例,如果換計算機把其中的1.1改乘其它即可,因為相信大家的cpu都應該在1G以上  ^_^
 cout<< (double) (timer/2.0/1000000000.0) << endl;