c54fft.cpp

VC下编程FFT 算法通过,可用在单片机,频谱分析效果良好.

#include "stdio.h"
#include "stdlib.h"
#include <Ctype.h>
#include <String.h>
#include <math.h>
//#define pi 3.1415926536
#define NN 128

void main(void)
 {int i,j,k,m,N,r,s,t;
 double pi= 4*atan(1);//3.1415926536;
 double a0,a1,b0,b1,p,z,g;
  int a[NN];
 double c[NN];
 double d[NN][2],w[NN][2];
 // char th1[15];
 //FILE *fp,*fp1;
 a[0]=0;
 a[1]=1;
 N=64;
 r=6;       //2**N

   for(i=0;i<N;i++)
   {c[i]=cos(2*1*pi*i/N);//+cos(2*2*pi*i/N);
   }
    for(i=0;i<N;i++)             //bit invers data
	{a[2*i]=a[i]/2;
      a[2*i+1]=a[2*i]+N/2;
	}
   g=0.0;
 //  fp=fopen("s1.dat","r");
   for(i=0;i<N;i++)    //  bit invers data
	{d[a[i]][0]=c[i];
      d[i][1]=0;
	  g=g+c[i];//求和
	}
 // fclose(fp);
 g=g/pow(2,r);//N   求平均值直流
    for(i=0;i<N;i++)
	{d[i][0]=d[i][0]-g;    //减直流
	}
    for(i=0;i<r;i++)
	{  
		for(m=0;m<pow(2,i);m++)//  generater sin 
		{w[m][0]=(cos(2*pi*m/pow(2,i+1)));
         w[m][1]=(-sin(2*pi*m/pow(2,i+1)));
		}
		for(j=0;j<N;j=(int)(j+pow(2,i+1)))
		{ 
			for(k=0;k<pow(2,i);k++)
			{
              s=(j+k);
			  t=(j+k+(int)pow(2,i));//碟形运算
			    a0=(d[s][0]+w[k][0]*d[t][0]-w[k][1]*d[t][1])/2;
				a1=(d[s][1]+w[k][0]*d[t][1]+w[k][1]*d[t][0])/2;
				b0=(d[s][0]-w[k][0]*d[t][0]+w[k][1]*d[t][1])/2;
				b1=(d[s][0]-w[k][0]*d[t][1]-w[k][1]*d[t][0])/2;
				d[s][0]=a0;    //存储real。
				d[s][1]=a1;    //imag
				d[t][0]=b0;     //real 
			  	d[t][1]=b1;     //imag
			}
		}
	}
for(k=0;k<N;k++)              //power 
   {p=(pow(d[k][0],2)+pow(d[k][1],2));
   z=sqrt(p);
  // a[k]=(int)1000*z;
   //z=d[k][0];
  // printf("%f\n",z);
   printf("%f\n",z);
 //  printf("%d\n",a[k]);
   }
}