FFT.zip_visual c

#include<stdio.h> #include<iostream> #include<stdlib.h> #include<math.h> #include"complex.h" const int N=2048; const int M=11; #define pi 3.141592654 #define e 2.71828 using namespace std; double max(double *array,int N) { double max_num=0; double *b=new double [N]; for(int i=0;i<N;i++) { b[i]=abs(array[i]); if(b[i]>max_num) max_num=b[i]; } return max_num; } int swap(int a,int M) { int number=0; for(int i=0;i<M;i++) { int temp=a%2; number=(number+temp)*2; a=a/2; } return number/2; } int find1(double *a,int N) { double max=0; int i; int position=0; int bound1=600.0/8000.0*N; int bound2=1000.0/8000.0*N; for(i=bound1;i<bound2;i++) { if(a[i]>max) { max=a[i]; position=i; } } return position; } int find2(double *a,int N) { double max=0; int i; int position=0; int bound1=1100.0/8000.0*N; int bound2=1700.0/8000.0*N; for(i=bound1;i<bound2;i++) { if(a[i]>max) { max=a[i]; position=i; } } return position; } void t2fft(int N,int M,complex array[]) { double P; for(int i=1;i<=M;i++) { int D=pow(double(2),(i-1)); for(int j=0;j<D;j++) { P=j*pow(double(2),M-i); double real=cos(2*pi*P/N); double imag=-sin(2*pi*P/N); complex rotate; rotate.set_complex(real,imag);//旋转因子； for(int K=j;K<N-1;K=K+pow(double(2),i)) { int index=(int)(K+D); complex temp=array[K]+array[index]*rotate; array[index]=array[K]-array[index]*rotate; array[K]=temp; } } } } int main() { int i; //用作循环计数 FILE *fp; fp=fopen("data1081.wav","rb");//为读，打开一个wav文件 if((fp=fopen("data1081.wav","rb"))==NULL) //若打开文件失败，退出 { printf("can't open this file\n"); exit(0); } fseek(fp,0,SEEK_SET); fseek(fp,0,SEEK_END); int f_size=ftell(fp);//文件长度； printf("此wav文件的长度为：%d\n",f_size); int data_length=(f_size-44)/2;//wav格式文件中有44个字节用来存储文件信息，故要滤除； printf("有效数据点数为：%d\n",data_length); //读取数据； double *data=new double [N]; fseek(fp,46,SEEK_SET); short temp;//wav文件存储数据类型为short； for(i=0;i<N;i++) { fread(&temp,sizeof(short),1,fp); data[i]=double(temp); } //处理数据过程； for(i=0;i<N;i++) { data[i]=(i>=data_length)?0:data[i];//做2048点FFT，不足点补零； } double max_num=max(data,N); for(i=0;i<N;i++) { data[i]=data[i]/max_num;//归一化； } //定义数组来存储FFT值； complex *FF=new complex [N]; double *F=new double [N];//存储频谱幅值； int note; for(i=0;i<N;i++) { note=swap(i,M); FF[i].set_complex(data[note],0);//按时间抽取，将数据按照倒序码排列； //F_imag[i]=0; } //FFT变换； t2fft(N,M,FF); for(i=0;i<N;i++) { //F[i]=sqrt(F_real[i]*F_real[i]+F_imag[i]*F_imag[i]); F[i]=FF[i].distance(); } //寻找两个峰值并输出； int x1; x1=find1(F,N); double f1=x1*8000/N; printf("第一个峰值的频率为：%f：\n",f1); int x2; x2=find2(F,N); double f2=x2*8000/N; printf("第二个峰值的频率为：%f：\n",f2); return 0; }