main.cpp

这是雨流技术法的vc++版本

#include"stdafx.h"
#include"rainflow.h"
#include"resource.h"

#include"mclmcr.h"
#include"mex.h"
#include<math.h>
#include<assert.h>

#define NMAX 5000000
#define H_TOLERANCE 0.005
#define PARIS_PARA 3.0

void push(double x[],int &ip,double t);
double get_x(double x[],int &index,int N);
void Four_Point_Judge(double temp[],int & ip,double y[],int &ny);

void rainflow(double x[],int nx,double y[],int & ny);

void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
	if(nrhs!=1)
		mexErrMsgTxt("one inputs requried.");
	if(nlhs>1)
		mexErrMsgTxt("Too many output arguments.");
	
	int nrows;
	int ncols;
	nrows=mxGetM(prhs[0]);
	ncols=mxGetN(prhs[0]);
	if(!mxIsDouble(prhs[0])||mxIsComplex(prhs[0])||nrows!=1)
		mexErrMsgTxt("Input must be a double vector");

	double *x;
	double *y;
	y=new double[NMAX];
	int ny=0;
	x=mxGetPr(prhs[0]);
//*********************************************
	rainflow(x,ncols,y,ny);//调用雨流计数函数
//*********************************************
	plhs[0]=mxCreateDoubleMatrix(1,ny,mxREAL);
	double *t;
	t=mxGetPr(plhs[0]);
	for(int i=0;i<ny;i++)
		t[i]=y[i];
	return;
}


void push(double x[],int &ip,double t)
{
	assert(ip!=NMAX&&ip>=0);
	x[ip]=t;
	ip++;
}

double get_x(double x[],int &index,int N)
{
	assert(index>=0&&index<N);
	index++;
	return x[index-1];
}

void Four_Point_Judge(double temp[],int & ip,double y[],int &ny)
{
	while(ip>=4)
	{
		if(temp[ip-1]>temp[ip-2]&&(temp[ip-1]-temp[ip-3])>=((-1)*H_TOLERANCE)&&(temp[ip-2]-temp[ip-4])>=((-1)*H_TOLERANCE))
		{
			push(y,ny,(temp[ip-3]-temp[ip-2])/2);
			temp[ip-3]=temp[ip-1];
			ip-=2;
			continue;
		}
		if(temp[ip-1]<temp[ip-2]&&(temp[ip-1]-temp[ip-3])<=H_TOLERANCE&&(temp[ip-2]-temp[ip-4])<=H_TOLERANCE)
		{
			push(y,ny,(temp[ip-2]-temp[ip-3])/2);
			temp[ip-3]=temp[ip-1];
			ip-=2;
			continue;
		}
		else
		{
			break;
		}
	
	}

	if(ip==3)
		if(fabs((temp[ip-1]-temp[ip-3]))<=H_TOLERANCE)
		{push(y,ny,fabs((temp[ip-1]-temp[ip-2])/2));ip-=2;
			return;}
}

void rainflow(double x[],int nx,double y[],int & ny)
{

	
	ny=0;
	int ip=0;
	//double temp[NMAX];
	double *temp;
	temp=new double[NMAX];
	double last,now,next;
	int index=0;
	last=get_x(x,index,nx);
	push(temp,ip,last);
	do{now=get_x(x,index,nx);}while(fabs(now-last)<H_TOLERANCE);//读取头两个不相等的应变
	while(index<nx)
	{
		next=get_x(x,index,nx);
		if(fabs(next-now)<H_TOLERANCE)
				continue;
		if(next>now&&now<last)
			{
				push(temp,ip,now);
				last=now;
				now=next;
				Four_Point_Judge(temp,ip,y,ny);
				continue;
			}
		if(next<now&&now>last)
			{
				push(temp,ip,now);
				last=now;
				now=next;
				Four_Point_Judge(temp,ip,y,ny);
				continue;
			}
		
		else
			{
				now=next;
				continue;
			}
	}
	if(index==nx)
		{
			push(temp,ip,now);
			Four_Point_Judge(temp,ip,y,ny);
		}
	int i=0;
	if(ip>=2)
	{	
		/*if(temp[0]<temp[1])
			i=1;
		for( ;i<ip;i+=2)       //第二次计数（峰值计数）
		{
			if(temp[i]>H_TOLERANCE)
				push(y,ny,temp[i]);
		}*/
		for(i=0;i<(ip-1);i++)       //第二次计数（峰值计数）
		{
			push(y,ny,fabs((temp[i+1]-temp[i])/2)/pow(2.0,1.0/PARIS_PARA));
		}
	}
			//push(y,ny,fabs((temp[i-1]-temp[i])/2)/pow(2.0,1.0/PARIS_PARA));

}