sound classify.cpp

神经网络中的声速测定的源代码 学神经网络的必用程序

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define train 120
#define test 120
#define PI (4*atan(1))

struct parameter
{
	int *data,emax;
	double energy,mean,delta,delta16,*mo,energy16,mo0;
};
	
struct netpara
{
	int iteration;
	double *doori,doorj,doork,doore,u1,u2,**wji,**wkj,c,door;
};

//定义存放图象信息的结构体内容
struct bmpinfo
{
 FILE *fp;
 long bfsize,bfoffbits,bisize,biwidth,biheight,bicompression,bisizeimage;
 int bibitcount,**blue,**red,**green;
};

void initial_para(struct parameter *para)
{
	para->data=NULL;
	para->mo=NULL;
}

//doori,doorj,doork,doore是输入层、隐层，输出层f()函数和误差总能量的门限值
//u1,u2是输出层对隐层和隐层对输入层的缩放系数
void initial_net(int N,struct netpara *para)
{
	int i;
	for(i=0;i<N;i++) para->doori[i]=0;
	para->wji=NULL;
	para->wkj=NULL;
	para->u1=0.05;
	para->u2=0.05;
	para->iteration=100;
	para->c=1;
	para->door=-1;
}

//释放一个2纬double型的矩阵
//n矩阵维数
//array矩阵指针
void freearray(int n,double **array)
{
	int i;
	for(i=0;i<n;i++) free(array[i]);
	free(array);
}

//释放颜色和灰度矩阵
void free_color_array(struct bmpinfo *bmpdata)
{
	int i;
	if(bmpdata->blue!=NULL)
	{
		for(i=0;i<bmpdata->biheight;i++)
		{
			free(bmpdata->blue[i]);
			free(bmpdata->green[i]);
			free(bmpdata->red[i]);
		}
		free(bmpdata->blue);
		free(bmpdata->green);
		free(bmpdata->red);
	}
	else
	{
		printf("\nFree_color_array: the color array has not been built!\n");
		exit(0);
	}
}

//以BMPINFO结构体里的信息为依据，写BMP文件头，无文件时创建，有同名的时覆盖
void writebmphead(struct bmpinfo *bmpdata,char outputbmp[20])
{
	FILE *out;
	long q=0;

	if((out=fopen(outputbmp,"wb+"))==NULL)
	{
		printf("\nsorry can not write head to the file %s!\n",outputbmp);
		exit(0);
	}
	
	rewind(out);
	putc(0x42,out);
	putc(0x4d,out);
	putw(bmpdata->bfsize,out);
	putw(0x00,out);
	putw(bmpdata->bfoffbits,out);
	putw(bmpdata->bisize,out);
	putw(bmpdata->biwidth,out);
	putw(bmpdata->biheight,out);
	putc(0x01,out);
	putc(0x00,out);
	putw(bmpdata->bibitcount,out);
	putw(bmpdata->bicompression,out);
	putw(bmpdata->bisizeimage,out);
	putw(q,out);
	putw(q,out);
	putw(q,out);
	putw(q,out);
	fclose(out);
}

//以BMPINFO结构体里的信息为依据，写BMP颜色信息，接着头文件写，不能创建文件
//写完后释放占用的空间
void writebmpcolor(struct bmpinfo *bmpdata,char outputbmp[20])
{
	long width;
	long i,j,k;
	FILE *out;
	if((out=fopen(outputbmp,"rb+"))==NULL)
	{
		printf("\nsorry can not write color to the file %s!\n",outputbmp);
		exit(0);
	}

	fseek(out,bmpdata->bfoffbits,0);
	
			width=bmpdata->biwidth*3;
			if (width/4*4!=width) width=(width/4+1)*4;
			for (i=bmpdata->biheight-1;i>=0;i--)
			{
				k=0;
				for(j=0;j<bmpdata->biwidth;j++)
				{
					k+=3;
					putc(bmpdata->blue[i][j],out);
					putc(bmpdata->green[i][j],out);
					putc(bmpdata->red[i][j],out);
				}
				for (;k<width;k++) putc(0x00,out);
			}
		//freefile(bmpdata);
		fclose(out);
			
}

//写BMP文件函数
void writebmp(struct bmpinfo *bmpdata,char outputbmp[20])
{
	writebmphead(bmpdata,outputbmp);
	writebmpcolor(bmpdata,outputbmp);
}

//门限函数f(x)
double f(double x)
{
	return(1/(1+exp(-x)));
}

//门限函数求导f'(x)
double df(double x)
{
	return(f(x)*(1-f(x)));
}

//读取原始数据，名称name必须含有统一通配符，加数字编号，从1开始
//datanum原始数据文件个数
//elementnum每个文件中提取数据个数
//name原始文件通配名称
//data原始数据存放的数组指针
//totalelement每个原始文件中实际提取数据个数
//int readdata(int elementnum,char name[20],int **data,int *totalelement)
void readdata(int datanum,int elementnum,char name[20],struct parameter *para,int *totalelement)
{
	FILE *ep;
	char inputfile[30];
	int i,j,k;
	struct parameter *add;
	
	add=para;
	for(i=0;i<datanum;i++)
	{
		sprintf(inputfile,"%s%d.dat",name,i+1);
		if((ep=fopen(inputfile,"rt+"))==NULL)
		{
			printf("\nSorry,the file name that you input: %s dose not exist!\n",inputfile);
			exit(0);
		}
		for(j=0;j<elementnum;j++)
			if(fscanf(ep,"%d",&add->data[j])==EOF) break;
		//将不够elementnum要求的数据补零
		if(j<elementnum)
			for(k=j;k<elementnum;k++) add->data[k]=0;

		totalelement[i]=j;

/*		printf("totalelement[%d]=%d\n",i,totalelement[i] );
		for(j=0;j<elementnum;j++)
			printf("data[%d][%d]=%d\n",i,j,data[i][j]);
*/		
		fclose(ep);
		add++;
	}
}	
//读取权值矩阵wji，wkj
void load_w(int N,int M,int L,struct netpara *para,char name[30])
{
	int i,j,k;
	FILE *ep;

	if((ep=fopen(name,"rt+"))==NULL)
	{
		printf("\nSorry can not read the file %s!\n",name);
		exit(0);
	}

	for(i=0;i<N;i++)
		if(fscanf(ep,"%f",&para->doori[i])==EOF) break;
	
	for(j=0;j<M;j++)
		for(i=0;i<(N+1);i++)
			if(fscanf(ep,"%f",&para->wji[j][i])==EOF) break;
	
	for(k=0;k<L;k++)
		for(j=0;j<(M+1);j++)
			if(fscanf(ep,"%f",&para->wkj[k][j])==EOF) break;
			
	fclose(ep);
}

//输出权值矩阵wji，wkj
void output_w(int N,int M,int L,struct netpara *para,char name[30])
{
	int i,j,k;
	FILE *ep;

	if((ep=fopen(name,"wt+"))==NULL)
	{
		printf("\nSorry can not write the file %s!\n",name);
		exit(0);
	}

	fprintf(ep,"\nThe door value of input level are:\n");

	for(i=0;i<N;i++) fprintf(ep,"\ndoori[%d]=%f",i,para->doori[i]);

	fprintf(ep,"\nThe best W corefference matrix is:\n");

	for(j=0;j<M;j++)
	{
		for(i=0;i<(N+1);i++)
			fprintf(ep,"\nwji[%d][%d]=%f",j,i,para->wji[j][i]);
		fprintf(ep,"\n");
	}
	
	for(k=0;k<L;k++)
	{
		for(j=0;j<(M+1);j++)
			fprintf(ep,"\nwkj[%d][%d]=%f",k,j,para->wkj[k][j]);
		fprintf(ep,"\n");
	}
	fclose(ep);
}

//正反fft函数，变换后的结果还存回原实虚矩阵
//N变换点数
//real实部矩阵
//image虚部矩阵
//if_fft判断是正变换还是反变换 1是正变换，非1是反变换
void fft(int N,double *real,double *image,int if_fft)
{
	int i,j,k,m,cache,no1,no2;
	double a,b,c=1,temp1r,temp1i,temp2r,temp2i;
	
	//码位倒读
	//得到碟形运算的级数，从1开始，比如N=8，那么m=3
	frexp(N,&m);
	m-=1;
	//采样点必须全部转换顺序
	for(i=0;i<N;i++)
	{
		k=i;
		cache=0;
		//按照手工10->2进制转换方法求模得到余数和结果，然后从低位到高位按照颠倒顺序的权值赋值
		for(j=0;j<m;j++)
		{
			cache+=(int)k%2*pow(2,m-1-j);
			k/=2;
		}
		//如果编号改变那么进行编号调换
		if(cache>i)
		{
			a=real[cache];
			b=image[cache];
			real[cache]=real[i];
			image[cache]=image[i];
			real[i]=a;
			image[i]=b;
		}
	}
	
	//对sin函数的符号进行赋值
	if(if_fft!=1) 
	{
//		printf("\nThe grade of the IFFT transform: %d\n",m);
		c=-1;
	}
	else 
	{
//		printf("\nThe grade of the FFT transform: %d\n",m);
	}

	//碟形运算
	//第几级
	//变换点数是N=2^m，所以有m级
	for(i=0;i<m;i++)
	{
		//第几群
		//第一级是N/2个群，之后每级减半，直到最后一级一个群
		for(j=0;j<pow(2,m-i-1);j++)
		{
			//第几个碟形变换
			//每级内需要计算的碟形变换数k从第一级的1个开始每级翻倍，直到最后一级的N/2个
			//每级内每次碟形运算将前级2个输入量经计算后转为本级对应编号的输入量
			//如果是IFFT的话按照x(n)=1/N*(FFT(X"(k)))"计算,"表示共轭运算符
			//其效果就是sin函数取负号
			for(k=0;k<pow(2,i);k++)
			{
				//得到本次碟形运算的2个输入量的编号，他们之间的差值从第一级的1开始逐级翻倍
				no1=pow(2,i+1)*j+k;
				no2=pow(2,i+1)*j+pow(2,i)+k;
				temp1r=real[no1];
				temp1i=image[no1];
				temp2r=real[no2];
				temp2i=image[no2];
				//W的方次递增步长由最后一级的1往前级逐级翻倍
				a=cos(2*PI/N*pow(2,m-i-1)*k);
				b=c*sin(2*PI/N*pow(2,m-i-1)*k);
 				real[no1]=temp1r+temp2r*a+temp2i*b;
				image[no1]=(temp1i-temp2r*b+temp2i*a);
				real[no2]=temp1r-temp2r*a-temp2i*b;
				image[no2]=(temp1i+temp2r*b-temp2i*a);
			}
		}
	}
	//如果是IFFT运算，那么所得结果x(n)前面要乘以系数1/N
	if(c==-1)
		for(i=0;i<N;i++)
		{
			real[i]/=N;
			image[i]/=N;
		}
}

//计算参数,只计算前N/2点的
//real,image是存放fft计算结果的数组
//所得结果存放在结构体para内
//para->mo存放的是每一位的模
//para->energy是所有模之和
//para->mean是所有位模的均值
//para->delta是所有位模的方差
//para->emax是所有位模值最大的点对应的坐标
void cal_parameter(int N,double *real,double *image,struct parameter *para)
{
	int i;
	double cache,mean;

	//求模和能量,忽略直流量
	if(para->mo==NULL) para->mo=(double*)malloc(N*sizeof(double));
	para->energy=0;
	for(i=1;i<N;i++)
	{
		para->mo[i]=pow(pow(real[i],2)+pow(image[i],2),0.5);
		para->energy+=para->mo[i];
//		printf("\npara->mo[%d]=%f",i,para->mo[i]);
	}
	para->energy16=0;
	for(i=1;i<17;i++) para->energy16+=para->mo[i];
	//均值
	para->mean=para->energy/(N-1);
	para->mo0=pow(pow(real[0],2)+pow(image[0],2),0.5);
	mean=para->energy16/16;
	//方差
	para->delta=0;
	for(i=1;i<N;i++) para->delta+=pow((para->mo[i]-para->mean),2)/(N-2);
	para->delta=pow(para->delta,0.5);

	para->delta16=0;
	for(i=1;i<17;i++) para->delta16+=pow((para->mo[i]-mean),2)/15;
	para->delta16=pow(para->delta16,0.5);

	//模值最大点的坐标
	cache=para->mo[1];
	para->emax=1;
	for(i=2;i<N-1;i++)
		if(cache<para->mo[i])
		{
			cache=para->mo[i];
			para->emax=i;
		}

	free(para->mo);
	
/*	printf("\npara->mo0=%f",para->mo0);
	printf("\npara->energy=%f",para->energy);
	printf("\npara->energy16=%f",para->energy16);
	printf("\npara->mean=%f",para->mean);
	printf("\npara->deltra=%f",para->delta);
	printf("\npara->deltra16=%f",para->delta16);
	printf("\npara->emax=%d",para->emax);
	printf("\nmax mo=%f",cache);
	printf("\n");
*/	
}

//BP神经网络
//samplenum是所有三类样本点个数
//fftnum是从每个样本文件中提取的数值个数,也是fft变换的点数
//iteration是循环迭代的最大次数
//N,M,L分别是输入层,隐层,输出层神经元个数
//para是存放所有样本fft变换之后所得数据的结构体数组指针
void bpnet(int samplenum,int fftnum,char name[20],int N,int M,int L,struct parameter *para,struct netpara *netpara)
{
	//judge是每类样本个数
	//totalelement是每个样本文件所得有用数据个数
	//class1，2，3是分辨结果每类样本个数
	int i,j,k,l,n,temp,judge,*totalelement,class1=0,class2=0,class3=0;
	//real，image是fft实部和虚部矩阵
	//netpj,netpk是隐层，输出层净输入
	//opi,opj,opk是输入层，隐层，输出层神经元的输出
	//awji,awkj是隐层和输出层对前级的权值调整量
	double *real,*image,**opi,**opj,**opk,**netpj,**netpk,**awkj,**awji;
	//ydk是期望输出量
	//enew,elast是总的误差能量
	//deltaj,deltak是隐层和输出层delta系数
	double *ydk,elast=0,enew=0,*deltaj,*deltak;
	struct parameter *add;
	double cache,c;
	struct bmpinfo bmpdata;
	
	bmpdata.bfoffbits=54;
	bmpdata.bisize=40;
	bmpdata.biheight=samplenum*5;
	bmpdata.biwidth=netpara->iteration/32;
	bmpdata.bicompression=0;
	bmpdata.bisizeimage=bmpdata.biwidth*bmpdata.biheight*3;
	bmpdata.bfsize=54+bmpdata.bisizeimage;
	bmpdata.bibitcount=24;

	//分配空间
	bmpdata.blue=(int **)malloc(bmpdata.biheight*sizeof(int *));
	bmpdata.green=(int **)malloc(bmpdata.biheight*sizeof(int *));
	bmpdata.red=(int **)malloc(bmpdata.biheight*sizeof(int *));
	for(i=0;i<bmpdata.biheight;i++)
	{
		bmpdata.blue[i]=(int *)malloc(bmpdata.biwidth*sizeof(int));
		bmpdata.green[i]=(int *)malloc(bmpdata.biwidth*sizeof(int));
		bmpdata.red[i]=(int *)malloc(bmpdata.biwidth*sizeof(int));
	}

	opi=(double**)malloc(samplenum*sizeof(double*));
	for(i=0;i<samplenum;i++)
		opi[i]=(double*)malloc((N+1)*sizeof(double));
	opj=(double**)malloc(samplenum*sizeof(double*));
	for(i=0;i<samplenum;i++)
		opj[i]=(double*)malloc((M+1)*sizeof(double));
	opk=(double**)malloc(samplenum*sizeof(double*));
	for(i=0;i<samplenum;i++)
		opk[i]=(double*)malloc(L*sizeof(double));

	ydk=(double*)malloc(L*sizeof(double));

	netpj=(double**)malloc(samplenum*sizeof(double*));
	for(i=0;i<samplenum;i++)
		netpj[i]=(double*)malloc(M*sizeof(double));
	netpk=(double**)malloc(samplenum*sizeof(double*));
	for(i=0;i<samplenum;i++)
		netpk[i]=(double*)malloc(L*sizeof(double));
	
	deltaj=(double*)malloc(M*sizeof(double));
	deltak=(double*)malloc(L*sizeof(double));
	
	awji=(double**)malloc(M*sizeof(double*));
	for(j=0;j<M;j++)
		awji[j]=(double*)malloc((N+1)*sizeof(double));
	
	awkj=(double**)malloc(L*sizeof(double*));
	for(j=0;j<L;j++)
		awkj[j]=(double*)malloc((M+1)*sizeof(double));
	
	//初始化
	for(i=0;i<bmpdata.biwidth;i++)
		for(j=0;j<bmpdata.biheight;j++)
		{
			
			if(j%30==0)
			{
				bmpdata.blue[j][i]=0;
				bmpdata.red[j][i]=0;
				bmpdata.green[j][i]=0;
			}
			else
				if(i%50==0)
				{
					bmpdata.blue[j][i]=0;
					bmpdata.red[j][i]=0;
					bmpdata.green[j][i]=0;
				}
				else
				{
					bmpdata.blue[j][i]=255;
					bmpdata.red[j][i]=255;
					bmpdata.green[j][i]=255;
				}
		}
	//权值，及权值调整量赋初值
	srand( (unsigned)time( NULL ) );
	for(j=0;j<M;j++)
		for(i=0;i<(N+1);i++)
		{
		//	srand( (unsigned)time( NULL ) );
			netpara->wji[j][i]=(double)rand()/32767;
		}
		
//	srand( (unsigned)time( NULL ) );
	for(i=0;i<L;i++)
		for(j=0;j<(M+1);j++)
		{
		//	srand( (unsigned)time( NULL ) );
			netpara->wkj[i][j]=(double)rand()/32767;
		}
		
/*	for(j=0;j<M;j++)
		for(i=0;i<(N+1);i++)
			printf("\nnetpara->wji[%d][%d]=%f",j,i,netpara->wji[j][i]);
	for(j=0;j<(M+1);j++)
		for(i=0;i<L;i++)
			printf("\nnetpara->wkj[%d][%d]=%f",i,j,netpara->wkj[i][j]);
*/
	//分配样本数据存储空间
	if(para->data==NULL)
	{
		add=para;
		for(i=0;i<samplenum;i++)
		{
			add->data=(int *)malloc(fftnum*sizeof(int));
			add++;
		}
	}
	
	totalelement=(int *)malloc(samplenum*sizeof(int));

	real=(double *)malloc(fftnum*sizeof(double));