caculateview.cpp

一套故障诊断监测系统

// CaculateView.cpp : implementation file
//

#include "stdafx.h"
#include "CDMS.h"
#include "CaculateView.h"
#include "math.h"//My code
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define PI 3.14159265
/////////////////////////////////////////////////////////////////////////////
// CCaculateView

IMPLEMENT_DYNCREATE(CCaculateView, CView)

CCaculateView::CCaculateView()
{
}

CCaculateView::~CCaculateView()
{
}


BEGIN_MESSAGE_MAP(CCaculateView, CView)
	//{{AFX_MSG_MAP(CCaculateView)
		// NOTE - the ClassWizard will add and remove mapping macros here.
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CCaculateView drawing

void CCaculateView::OnDraw(CDC* pDC)
{
	CDocument* pDoc = GetDocument();
	// TODO: add draw code here
}

/////////////////////////////////////////////////////////////////////////////
// CCaculateView diagnostics

#ifdef _DEBUG
void CCaculateView::AssertValid() const
{
	CView::AssertValid();
}

void CCaculateView::Dump(CDumpContext& dc) const
{
	CView::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CCaculateView message handlers

//获得数组的最大值,buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetMax(double buffer[], int length)
{
	int i;
	double temp=buffer[0];
	for(i=1;i<length;i++)
	{
		if(temp<buffer[i])
			temp=buffer[i];
	}
	return temp;
}

//获得数组的最小值,buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetMin(double buffer[], int length)
{
	int i;
	double temp=buffer[0];
	for(i=1;i<length;i++)
	{
		if(temp>buffer[i])
			temp=buffer[i];
	}
	return temp;
}
//获得数组的峰峰值,buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetPeakPeakValue(double buffer[], int length)
{
	double max, min, temp;

	max = GetMax(buffer, length);
	min = GetMin(buffer, length);
	temp = max - min;	

	return temp;
}
//获得数组的平均值,buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetAverage(double buffer[], int length)
{
	double sum=0;
	for(int i=0;i<length;i++)
	{
		sum=sum+buffer[i];
	}
	double average=sum/(double)length;
	return average;
}
//获得数组的均方值，buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetAverSquare(double buffer[], int length)
{
	double sum=0;
	for(int i=0;i<length;i++)
	{
		sum=sum+buffer[i]*buffer[i];
	}
	double average=sum/(double)length;
	return average;
}
//获得数组的有效值（均方根值）,buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetVirtualValue(double buffer[], int length)
{
	double temp=0;
	for(int j=0;j<length;j++)
	{
		temp=temp+buffer[j]*buffer[j];			
	}
	temp=sqrt(temp/(double)length);

	return temp;
}
//////////////////////////////////////////////////
//获得数组的均方幅值，均方幅值的大小为均方值的开方
//////////////////////////////////////////////////
//获得数组的方差，Average为平均值，AverSquare为均方值
double CCaculateView::GetSquDiffer(double Average, double AverSquare)
{
	double m_dSquDiffer=0;
	m_dSquDiffer=AverSquare-Average*Average;
	return m_dSquDiffer;
}

//获得数组的峭度，buffer[]为数组的指针，length为数组的长度,AverSquare为均方值
double CCaculateView::GetKurtosis(double buffer[], int length,double AverSquare)
{
	double sum=0;
	for(int i=0;i<length;i++)
	{
		sum=sum+pow(buffer[i],4);
	}
	sum=sum/(double)length;
	sum=sum/(AverSquare*AverSquare);
	return sum;
}
//获得数组的平均幅值，buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetAverAm(double buffer[], int length)
{
	double sum=0;
	for(int i=0;i<length;i++)
	{
		sum=sum+fabs(buffer[i]);
	}
	sum=sum/length;
	return sum;
}
//获得数组的方根幅值，buffer[]为数组的指针，length为数组的长度
double CCaculateView::GetExtractAm(double buffer[],int length)
{
	double sum=0;
	for(int i=0;i<length;i++)
	{
		sum=sum+sqrt(fabs(buffer[i]));
	}
	sum=sum/length;
	sum=sum*sum;
	return sum;
}
////////////////////////////////////////////////////////////////
//////峰值=最大值和最小值中绝对较大值
//////波形指标=均方幅值/平均幅值
//////峰值指标=峰值/均方幅值
//////脉冲指标=峰值/平均幅值
//////裕度指标=峰值/方根幅值
////////////////////////////////////////////////////////////////
/*
关于FFT变换函数参数的说明
  pr——双精度实型一维数组，长度为n。
		当l=0时，存放n个采样输入的实部，返回时存放离散傅立叶变换的模；
        当l=1时，存放傅立叶变换的n个实部，返回时存放逆傅立叶变换的模。
  pi——双精度实型一维数组，长度为n。
		当l=0时，存放n个采样输入的虚部，返回时存放离散傅立叶变换的幅角；
        当l=1时，存放傅立叶变换的n个虚部，返回时存放逆傅立叶变换的幅角。
		其中幅角的单位为度。
  n——整型变量。输入的点数。
  k——整型变量。满足n=2（K）。
  fr——双精度实型一维数组，长度为n。
		当l=0时，返回傅立叶变换的实部；
        当l=1时，返回逆傅立叶变换的实部。
  fi——双精度实型一维数组，长度为n。
		当l=0时，返回傅立叶变换的虚部；
        当l=1时，返回逆傅立叶变换的虚部。
  l——整型变量。
		若l=0，表示要求本函数计算傅立叶变换。
       若l=1，表示要求本函数计算逆傅立叶变换。
  il——整型变量。
		若il=0，表示不要求本函数计算傅立叶变换或逆傅立叶变换的模与幅角。
        若il=1，表示要求本函数计算傅立叶变换或逆傅立叶变换的模与幅角。
*/
void CCaculateView::FFT(double pr[], double pi[], int n, int k,
						double fr[], double fi[], int l, int il)
{
	int it,m,is,i,j,nv,l0;
	double p,q,s,vr,vi,poddr,poddi;
	for(it=0;it<=n-1;it++)
	{
		m=it;
		is=0;
		for(i=0;i<=k-1;i++)
		{
			j=m/2;
			is=2*is+(m-2*j);
			m=j;
		}
		fr[it]=pr[is];
		fi[it]=pi[is];
	}
	pr[0]=1.0;
	pi[0]=0.0;
	p=6.283185306/(1.0*n);
	pr[1]=cos(p);
	pi[1]=-sin(p);
	if(l!=0)
		pi[1]=-pi[1];
	for(i=2;i<=n-1;i++)
	{
		p=pr[i-1]*pr[1];
		q=pi[i-1]*pi[1];
		s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
		pr[i]=p-q;
		pi[i]=s-p-q;
	}
	for(it=0;it<=n-2;it=it+2)
	{
		vr=fr[it];
		vi=fi[it];
		fr[it]=vr+fr[it+1];
		fi[it]=vi+fi[it+1];
		fr[it+1]=vr-fr[it+1];
		fi[it+1]=vi-fi[it+1];
	}
	m=n/2;
	nv=2;
	for(l0=k-2;l0>=0;l0--)
	{
		m=m/2;
		nv=2*nv;
		for(it=0;it<=(m-1)*nv;it=it+nv)
			for(j=0;j<=(nv/2)-1;j++)
			{
				p=pr[m*j]*fr[it+j+nv/2];
				q=pi[m*j]*fi[it+j+nv/2];
				s=pr[m*j]+pi[m*j];
				s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
				poddr=p-q;
				poddi=s-p-q;
				fr[it+j+nv/2]=fr[it+j]-poddr;
				fi[it+j+nv/2]=fi[it+j]-poddi;
				fr[it+j]=fr[it+j]+poddr;
				fi[it+j]=fi[it+j]+poddi;
			}
	}
	if(l!=0)
	{
		for(i=0;i<=n-1;i++)
		{
			fr[i]=fr[i]/(1.0*n);
			fi[i]=fi[i]/(1.0*n);
		}
	}

	for(i=0;i<n;i++)
	{	
		fr[i]=fr[i]/n;
		fi[i]=fi[i]/n;
		if(il!=0)
		{
			pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i])*2;
			pi[i]=atan2(fi[i],fr[i])*180.0/3.141592653589793;
			if(pi[i]<0)
				pi[i]+=360.0;
		}
	}
	return;
}
//设fRPM为转子转速，ab为数据指针，k为特征频谱倍数。
//设fSample为采样频率。
//DFT算法：(num=60*fSample/fRPM;当求转频幅值，k=1；当求2倍频幅值