hht.cpp

希尔伯特—黄变换的控制台程序

#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>
#include<math.h>
#include<iostream.h>
/*************************************************************
读数据，将硬盘数据读入内存。
**************************************************************/
void ReadDat(int array_length,float y[])
{
	int i;
	FILE *fp;
	if((fp=fopen("newin.dat","r"))==NULL)
	{
		cout<<"can not open file!\n";
		return;
	}
	for(i=0;i<array_length;i++)
		if(!feof(fp))
			fscanf(fp,"%f,",&y[i]);
		else
		{
			cout<<"end of data!";
			fclose(fp);
			return;
		}
}
/*****************************************************************
extract_max()为求极大值函数,y[]为被求数据，array_max[]中存放结果。
******************************************************************/
void extract_max(double y[],int array_length,double array_max[])
{  
   int i=0,j=0;
   for(j=0,i=0;i<array_length-2;i++)
       if((y[i+1]>y[i])&&(y[i+1]>y[i+2]))
       {
          array_max[j+1]=y[i+1];
          j++;
       }
	   else
	   {
		  array_max[j+1]=0;
		  j++;
	   }
	array_max[0]=0;
	array_max[array_length-1]=0;
}
/****************************************************************
extract_min()为求极小值函数,同上。
*****************************************************************/
void extract_min(double y[],int array_length,double array_min[])
{  
   int i=0,j=0;
   for(j=0,i=0;i<array_length-2;i++)
       if((y[i+1]<y[i])&&(y[i+1]<y[i+2]))
       {
          array_min[j+1]=y[i+1];
          j++;
       }
       else
       {
          array_min[j+1]=0;
          j++;
       }
	array_min[0]=0;
	array_min[array_length-1]=0;
}
/*******************************************************
nnz()计算一维数组中非零元个数。
*********************************************************/
int nnz(double y[],int array_length)
{
	int i=0,nnz=0;
	for(i=0;i<array_length;i++)
		if(y[i]!=0)nnz++;
	return nnz;
}
/********************************************************
找出数据y[]中非零元下标,存入数组x[]。
*********************************************************
void index_notzero( double y[],int array_length,int x[])
{
	int i=0,j=0;
	for(i=0;i<array_length;i++)
		if(y[i]!=0)
		{
			x[j]=i;
			j++;
		}
}
/*****************************************************
三次样条插值核心函数。
入口参数：
        x：给定的节点；
		y：给定节点处的对应的值；
		n：已有节点个数，n≥3；
		xx：插值点；
		m：插值点数；
出口参数：
        yy：插值结果；
******************************************************/
void SP(double*x,double*y,int n,double*xx,int m,double*yy)
{
	int i,j,k;
	double z,h1,h2,h3,h4;
	double *h,*d,*s,*s2,*e;
	h=(double*)calloc(n,sizeof(double));
	if(h==NULL)
		exit(1);
	d=(double*)calloc(n,sizeof(double));
	if(d==NULL)
		exit(1);
	s=(double*)calloc(n,sizeof(double));
	if(s==NULL)
		exit(1);
	s2=(double*)calloc(n,sizeof(double));
	if(s2==NULL)
		exit(1);
	e=(double*)calloc(n,sizeof(double));
	if(e==NULL)
		exit(1);
	for(i=0;i<=n-2;i++)
	{
		h[i]=x[i+1]-x[i];
		d[i]=(y[i+1]-y[i])/h[i];
	}
	s2[0]=s2[n-1]=0.0;
	for(i=1;i<=n-2;i++)
    s2[i]=6.0*(d[i]-d[i-1]);
	z=0.5/(h[0]+h[1]);
	s[0]=-h[1]*z;
	e[0]=s2[1]*z;
	for(i=1;i<=n-3;i++)
	{
		z=1.0/(2.0*(h[i]+h[i+1])+h[i]*s[i-1]);
		s[i]=-h[i+1]*z;
		e[i]=(s2[i+1]-h[i]*e[i-1])*z;
	}
	s2[n-2]=e[n-3];
	for(i=n-3;i>=1;i--)
		s2[i]=s[i-1]*s2[i+1]+e[i-1];
	for(i=0;i<=n-2;i++)
		s[i]=(s2[i+1]-s2[i])/h[i];
	i=1;
	k=0;
	for(j=0;j<=m-1;j++)
	{
		if(xx[j]>x[i])
		{
			k=i;
			i+=1;
		}
	h1=xx[j]-x[k];
	h2=xx[j]-x[i];
	h3=h1*h2;
	h4=s2[k]+h1*s[k];
	z=(s2[i]+s2[k]+h4)/6.0;
	yy[j]=y[k]+h1*d[k]+h3*z;
	}
	free(h);
	free(d);
	free(s);
	free(s2);
	free(e);
}
/**************************************************************************
三次线性内插。
***************************************************************************/
void spline(double array[],int array_length,double yy[])//含零的array[]!
{
	int i=0,j=0,n=0;
	double *x;
	double *y;
	int m;
	double *xx;
	double time_begin=0.0,time_end=10.0;
	
	double step;
	/*************************************
    SP()入口参数：极值点数目。
    **************************************/
	n=nnz(array,array_length);
    /*************************************
    SP()入口参数：*x,极值点对应的横坐标。
    **************************************/
	x=(double*)calloc(n,sizeof(double));
	if(x==NULL)
		exit(1);
	for(j=0,i=0;i<array_length;i++)
		if(array[i]!=0)
		{
			x[j]=((i+1)/((double)array_length))*(time_end-time_begin);
			//cout<<x[j]<<',';
			j++;
			
		}
    /*************************************
    SP()入口参数：极值点：y[]。
    //存放极值的数组（除去上述为零处）
    **************************************/
	y=(double*)calloc(n,sizeof(double));
	if(y==NULL)
		exit(1);
	for(j=0,i=0;i<array_length;i++)
		if(array[i]!=0)
		{
			y[j]=array[i];
			j++;
		}
		//cout<<j<<endl;
    /*******************************
    SP()入口参数：要插值点的数目。
    ********************************/
	m=array_length;
    /**********************************
    SP()入口参数：xx[]准备插值处(横坐标)。
    ***********************************/
	xx=(double*)calloc(m,sizeof(double));
	if(xx==NULL)
		exit(1);

	step=(x[n-1]-x[0])/((double)m-1.0);
	xx[0]=x[0];
	xx[m-1]=x[n-1];
	for(i=1;i<m;i++)
	{
		if((xx[i-1]+step)<=x[n-1])
			xx[i]=xx[i-1]+step;
		//cout<<xx[i]<<',';
	}
	//cout<<xx[511]<<endl;
	//cout<<"xx[511]="<<xx[511]<<endl;ok!
	/**********插值********/
	SP(x,y,n,xx,m,yy);
	free(x);
	free(y);
	free(xx);

}
/********************************************************************************
提取固有模态函数。
**********************************************************************************/
void extract_imf(double y[],double yy_max[],double yy_min[],int array_length,double h[])
{
	int i=0;
	for(i=0;i<array_length;i++)
	{
		h[i]=y[i]-(yy_max[i]+yy_min[i])/2.0;
		//cout<<h[i]<<'h';ok!
	}
}
/********************************************************************************
判断是否是固有模态函数。
**********************************************************************************/
int isimf(double h[],int array_length)
{
	int i=0,k1=0,k2=0,zeros=0,extremums=0;
	double h_max=0,m_max=0;
    double *array_max,*array_min;
	int n_max=0,n_min=0;
	double *yy_max,*yy_min;
	double *m;
	/******首先计算零点*************/
	for(zeros=0,i=0;i<array_length-1;i++)
		if(h[i]*h[i+1]<=0)
			zeros+=1;
	/******************近似算法*************/
		for(i=0;i<array_length;i++)
		{
			if(fabs(h[i])>h_max)
				h_max=fabs(h[i]);
		}
	/**固有模态函数的第一个条件**/
	/**重复main()中求包络的过程*/
    array_max=(double*)calloc(array_length,sizeof(double));
	if(array_max==NULL)
		exit(1);
	array_min=(double*)calloc(array_length,sizeof(double));
	if(array_min==NULL)
		exit(1);
	extract_max(h,array_length,array_max);
    extract_min(h,array_length,array_min);
	/***********************************************
	FILE *fp;
	fp=fopen("array_max.dat","w");
	for(i=0;i<array_length;i++)
	{
		fprintf(fp,"%f%c",array_max[i],',');
	}
	fp=fopen("array_min.dat","w");
	for(i=0;i<array_length;i++)
	{
		fprintf(fp,"%f%c",array_min[i],',');
	}
	fclose(fp);
	************************************************/

	n_max=nnz(array_max,array_length);
	n_min=nnz(array_min,array_length);

	yy_max=(double*)calloc(array_length,sizeof(double));
	if(yy_max==NULL)
		exit(1);
	yy_min=(double*)calloc(array_length,sizeof(double));
	if(yy_min==NULL)
		exit(1);

	spline(array_max,array_length,yy_max);
	spline(array_min,array_length,yy_min);
	//for(i=0;i<array_length;i++)
	//	cout<<yy_max[i]<<',';

    /*******************test**********************
	fp=fopen("yy_max.dat","w");
	for(i=0;i<array_length;i++)
	{
		fprintf(fp,"%f%c",yy_max[i],',');
	}
	fp=fopen("yy_min.dat","w");
	for(i=0;i<array_length;i++)
	{
		fprintf(fp,"%f%c",yy_min[i],',');
	}
	fclose(fp);
    /******************test**********************/
	free(array_max);
	free(array_min);
    
	m=(double*)calloc(array_length,sizeof(double));
	if(m==NULL)
		exit(1);
	for(i=0;i<array_length;i++)
	{
		m[i]=fabs((yy_max[i]+yy_min[i])/2);
	}
	for(i=0;i<array_length;i++)
	{
		if(fabs(m[i])>m_max)
			m_max=fabs(m[i]);
	}
	if(m_max<h_max/10000)
	{
		k1=1;
		cout<<"条件一已满足!"<<endl;
	}
	else
		k1=0;
	
	/**固有模态函数的第二个条件**/
	
	extremums=n_max+n_min;

	//cout<<"极点:"<<extremums<<endl;
	//cout<<"零点:"<<zeros<<endl;

	if(abs(extremums-zeros)<=1)
		k2=1;
	else
		k2=0;
	
	
	free(yy_max);
	free(yy_min);
	free(m);
	//for(i=0;i<array_length;i++)
	//	cout<<h[i]<<',';
	//cout<<endl;
	return (k1&&k2);
}
/***************************************************************************************
主函数************************************************************************************
******************************************************************************************
**************************************************************************************/
void main()
{
	int i=0,j=0;
	double *array_max,*array_min;
	int n_max=0,n_min=0;
	double *yy_max,*yy_min;
	double *h;//*h1;
	int isorno=0;
	int test=0;
	double *y;
	float *y1;
	int array_length=512;
	y=(double*)calloc(array_length,sizeof(double));
	if(y==NULL)
		exit(1);
    y1=(float*)calloc(array_length,sizeof(float));
	if(y==NULL)
		exit(1);
	ReadDat(array_length,y1);
	for(i=0;i<array_length;i++)
		cout<<y1[i]<<',';
	for(i=0;i<array_length;i++)
		y[i]=(double)y1[i];
	free(y1);
	/***********************************
     计算原始一维地震数据长度
    ************************************
    array_length=sizeof(y)/sizeof(y[0]);
    /***********************************
	参数初始化。
    ***********************************/
notimf:
	n_max=0,n_min=0;
	isorno=0;
	//存放极值的数组（非极值处为零）
	array_max=(double*)calloc(array_length,sizeof(double));
	if(array_max==NULL)
		exit(1);
	array_min=(double*)calloc(array_length,sizeof(double));
	if(array_min==NULL)
		exit(1);
	/*****************************
	提取极值，但非极值处用0填充。
	******************************/
	extract_max(y,array_length,array_max);
    extract_min(y,array_length,array_min);
	/*************************************
    极值点数目：n_max,n_min。
    **************************************/
	//n_max=nnz(array_max,array_length);
	//n_min=nnz(array_min,array_length);
	//cout<<n_max<<endl<<n_min;ok!
    /*********************************************
    初始化SP()出口参数：yy_max[],yy_min[]，存放插值结果。
    **********************************************/
	yy_max=(double*)calloc(array_length,sizeof(double));
	if(yy_max==NULL)
		exit(1);
	yy_min=(double*)calloc(array_length,sizeof(double));
	if(yy_min==NULL)
		exit(1);
	spline(array_max,array_length,yy_max);
	spline(array_min,array_length,yy_min);
	free(array_max);
	free(array_min);
	/*double max=0,min=0;
	for(i=0;i<array_length;i++)
		if(max<fabs(yy_max[i]))
			max=yy_max[i];
			cout<<"yy_max="<<max<<','<<endl;
	for(i=0;i<array_length;i++)
		if(min>(yy_max[i]))
			min=yy_max[i];
			cout<<"yy_min="<<min<<','<<endl;
			cout<<"426="<<yy_min[426]<<endl;
    /***********************************
    以下程序将插值结果的数据输出到文件。
    ************************************
	FILE *fp;
	fp=fopen("yy_max.dat","w");
	for(i=0;i<array_length;i++)
	{
		fprintf(fp,"%f%c",yy_max[i],',');
	}
	fp=fopen("yy_min.dat","w");
	for(i=0;i<array_length;i++)
	{
		fprintf(fp,"%f%c",yy_min[i],',');
	}
	fclose(fp);
	/***********************************************
    提取固有模态函数，数据放入数组h[]。
    ************************************************/
	h=(double*)calloc(array_length,sizeof(double));
	if(h==NULL)
		exit(1);
	//h1=(double*)calloc(array_length,sizeof(double));
	//if(h1==NULL)//因为h[i]的值在函数调用时被改变!
	//	exit(1);
	extract_imf(y,yy_max,yy_min,array_length,h);
	free(yy_max);
	free(yy_min);
	/*for(i=0;i<array_length;i++)
		cout<<h[i]<<',';
	cout<<endl;

	/*
	double max=0;
	for(i=0;i<array_length;i++)
		if(max<fabs(h[i]))
			max=h[i];
			cout<<"h_max="<<max<<','<<endl;
	*/
	/*******************************************************/
	isorno=isimf(h,array_length);
	if(isorno)
	{
		cout<<"找到固有模态函数!并存入h1.dat中!"<<endl;
		cout<<"求了"<<test<<"次!"<<endl;
	}
	if(isorno==0)
	{
		cout<<"不是固有模态函数!"<<endl;
    	test+=1;
		for(i=0;i<array_length;i++)
		{

			y[i]=h[i];
		}
		//for(i=0;i<array_length;i++)
		//cout<<h[i]<<','<endl;
	}
	else
	{
		FILE *fp1;
		fp1=fopen("h1.dat","w");
		for(i=0;i<array_length;i++)
		{
			fprintf(fp1,"%f%c",h[i],',');
		}
		fprintf(fp1,"%c",'\n');
		fclose(fp1);
		free(h);
		return;
	}
	if(test>=999999)
	{
	    FILE *fp1;
		fp1=fopen("h1.dat","w");
		for(i=0;i<array_length;i++)
		{
			fprintf(fp1,"%f%c",h[i],',');
		}
		fprintf(fp1,"%c",'\n');
		fclose(fp1);
		free(h);
		return;
	}
	free(h);
	goto notimf;
}