#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;
}
