signalcalc.cpp
来自「常用算法与数据结构原代码」· C++ 代码 · 共 1,624 行 · 第 1/3 页
CPP
1,624 行
// SignalCalc.cpp : Defines the initialization routines for the DLL.
//
#include "stdafx.h"
#include "SignalCalc.h"
#include <math.h>
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#define PI 3.14159265
#define max_nscales 20
#define max_wlen 23
#define NODS 20
//
// Note!
//
// If this DLL is dynamically linked against the MFC
// DLLs, any functions exported from this DLL which
// call into MFC must have the AFX_MANAGE_STATE macro
// added at the very beginning of the function.
//
// For example:
//
// extern "C" BOOL PASCAL EXPORT ExportedFunction()
// {
// AFX_MANAGE_STATE(AfxGetStaticModuleState());
// // normal function body here
// }
//
// It is very important that this macro appear in each
// function, prior to any calls into MFC. This means that
// it must appear as the first statement within the
// function, even before any object variable declarations
// as their constructors may generate calls into the MFC
// DLL.
//
// Please see MFC Technical Notes 33 and 58 for additional
// details.
//
/////////////////////////////////////////////////////////////////////////////
// CSignalCalcApp
BEGIN_MESSAGE_MAP(CSignalCalcApp, CWinApp)
//{{AFX_MSG_MAP(CSignalCalcApp)
// NOTE - the ClassWizard will add and remove mapping macros here.
// DO NOT EDIT what you see in these blocks of generated code!
//}}AFX_MSG_MAP
END_MESSAGE_MAP()
/////////////////////////////////////////////////////////////////////////////
// CSignalCalcApp construction
CSignalCalcApp::CSignalCalcApp()
{
// TODO: add construction code here,
// Place all significant initialization in InitInstance
}
/////////////////////////////////////////////////////////////////////////////
// The one and only CSignalCalcApp object
CSignalCalcApp theApp;
int CalculatePower(int N)
{
int i,m;
for(i=0;i<20;i++)
{
m=i+1;
N/=2;
if(N==1) break;
}
return m;
}
/* NOTE: WINDOWFLAG=1 no window
WINDOWFLAG=2 hanning window
WINDOWFLAG=3 hamming window
WINDOWFLAG=4 blackman window */
void windowed( float *x, int n, int windowflag )
{ float multiplier;
int i;
switch(windowflag)
{ case 1:break;
case 2:
for(i=0;i<n;i++)
{ multiplier=(float)(0.5*(1-cos(6.2831853*i/(n-1))));
x[i]=x[i]*multiplier;
}
break;
case 3:
for(i=0;i<n;i++)
{ multiplier=(float)(0.54-0.46*cos(6.2831853*i/(n-1)));
x[i]=x[i]*multiplier;
}
break;
case 4:
for(i=0;i<n;i++)
{ multiplier=(float)(0.42-0.5*cos(2*PI*i/(n-1))+0.08*cos(4*PI*i/(n-1)));
x[i]=x[i]*multiplier;
}
break;
}
}
/////////////////峭度系数//////////////////////////////////////////////
extern "C" __declspec(dllexport)float kurtosis(float *x, int N)// 采用MATLAB中的峭度计算公式 smj
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
float ans;
float aver = 0.,sum1=0.,sum2=0.;
for(int i=0;i<N;i++)
aver += x[i];
aver = aver/N;//样本序列均值aver
for(i=0;i<N;i++)
{
ans=x[i]-aver;
sum1+=ans*ans*ans*ans;
sum2+=ans*ans;
}
return sum1*N/(sum2*sum2);//返回峭度系数值
}
///////////////////////////////////////////////////////////////////////
// 快速傅立叶变换
// xr--输入序列的实部；xi--输入序列的虚部；
// m--序列长度N取2的对数，即 ；
// inv--0为正变换，1为反变换；
extern "C" __declspec(dllexport) void fft(float *xr,float *xi,int m,int inv)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
float ur,ui,wr,wi,tr,ti,tmr,tmi;
int n,l,i,j,le,le1,ip,nv1,nv2,k;
n=(int)pow(2.0,(double)m);
//calculate fft
for(l=1;l<=m;l++)
{
le=(int)pow(2.0,(m+1-l));
le1=le/2;
ur=1.0f;
ui=0.0f;
wr=(float)cos(PI/le1);
if(inv==0)
wi=-(float)sin(PI/le1);
else
wi=(float)sin(PI/le1);
for(j=0;j<le1;j++)
{
for(i=j;i<n;i+=le)
{
ip=i+le1;
tmr=xr[i]-xr[ip];
tmi=xi[i]-xi[ip];
xr[i]=xr[i]+xr[ip];
xi[i]=xi[i]+xi[ip];
xr[ip]=tmr*ur-tmi*ui;
xi[ip]=tmr*ui+tmi*ur;
}
if((j+1)==le/4) //to enhance precision
{
if(inv==0)
{
ur=0.0f;
ui=-1.0f;
}
else
{
ur=0.0f;
ui=1.0f;
}
}
else
{
tr=ur*wr-ui*wi;
ui=ur*wi+ui*wr;
ur=tr;
}
}
}
// reverse order
nv2=n/2;
nv1=n-1;
j=0;
for(i=0;i<nv1;i++)
{ if(i<j)
{ tr=xr[j];
ti=xi[j];
xr[j]=xr[i];
xi[j]=xi[i];
xr[i]=tr;
xi[i]=ti;
}
k=nv2;
while(k<=j)
{ j-=k;
k=k/2;
}
j+=k;
}
if(inv==1)
{ for(i=0;i<=n-1;i++)
{ xr[i]=xr[i]/n;
xi[i]=xi[i]/n;
}
}
}
/////////////////////////////////////////////
void fht2t( float *x, int n, int m )
{ int i,j,k,n1,n2,n4,l1,l2,l3,l4;
float xt,ophase1,ophase2,ss,cc,w,a;
// digit reverse counter,use radar algorithm
j=0;
for(i=0;i<n-1;i++)
{ if(i<j)
{ xt=x[i]; x[i]=x[j]; x[j]=xt; }
k=n>>1; //k=n/2
while(k<=j)
{ j-=k;
k=k>>1; //k=k/2
}
j=j+k;
}
//Hartley transform
for( i=0;i<n;i+=2)
{ xt=x[i];
x[i]=xt+x[i+1];
x[i+1]=xt-x[i+1];
}
n2=1;
for( k=1;k<m;k++ )
{ n4=n2;
n2=n4<<1; //n2=n4*2
n1=n2<<1; //n1=n2*2
w=(float)(6.2831853/n1);
for( j=0;j<n;j+=n1 )
{ l2=j+n2;
l3=j+n4;
l4=l2+n4; //j=0;l2=n/2,l3=n/4,l4=3*n/4
xt=x[j];
x[j]=xt+x[l2];
x[l2]=xt-x[l2];
xt=x[l3];
x[l3]=xt+x[l4];
x[l4]=xt-x[l4];
a=w;
for( i=1; i<n4;i++ )
{ l1=j+i;
l2=l1+n2;
l3=j+n2-i;
l4=l3+n2; //l1=k,l2=n/2+k,l3=n/2-k,l4=n-k
ss=(float)sin(a);
cc=(float)cos(a);
ophase1=x[l2]*cc+x[l4]*ss;
ophase2=x[l2]*ss-x[l4]*cc;
xt=x[l1];
x[l1]=xt+ophase1;
x[l2]=xt-ophase1;
xt=x[l3];
x[l3]=xt+ophase2;
x[l4]=xt-ophase2;
a+=w;
}
}
}
}
////////////////////////////////////////////////////////////////
/* NOTE: outport order: */
/* {Re[0],Re[1],...Re[N/2],Im[N/2-1],Im[N/2-2]...Im[1] */
/////////////////////////////////////////////////////////////////
extern "C" __declspec(dllexport)void efft( float *x, int n, int windowflag )
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
float xt;
int i,m;
m=CalculatePower(n);
windowed( x, n, windowflag);
fht2t( x, n, m );
for(i=1;i<n/2;i++)
{ xt=x[i];
x[i]=(xt+x[n-i])/2;
x[n-i]=-(xt-x[n-i])/2;
}
}
//////////////////////校正的频谱计算/////////////////////////////////////////
// x--输入实序列，输出的前N/2+1项为校正后的幅值谱（bRectify=FALSE则不校正）；
// N--序列长度；
// samplefrequence--采样频率；
// windowflag--加窗类型，1为矩形窗，2为汉宁窗，3为海明窗，4为布莱克曼窗；
// spectrumpeaknumber--需校正的幅值谱数目；
extern "C" __declspec(dllexport) void amendfft( float *x,int N,float samplefrequence,
int windowflag,BOOL bRectify,int spectrumpeaknumber)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,j,k = 0, l;
float max,amplititude,phase,frequence,kt,aver=0.0f;
float *save, *y;
HANDLE hFloatsave,hFloaty;
float ophase1,ophase2,dk;
max=0.0f;
hFloaty=(HANDLE)GlobalAlloc(GMEM_MOVEABLE|GMEM_ZEROINIT,N*sizeof(float));
y=(float *)GlobalLock(hFloaty);
//原始数据去均值2001/03/06
for(i=0;i<N;i++)
aver+=x[i];
aver=aver/N;
for(i=0;i<N;i++)
x[i]-=aver;
if(windowflag>2)
{ hFloatsave=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
save=(float *)GlobalLock(hFloatsave);
for(i=0;i<N/2;i++)
save[i]=x[i];
for(i=N/2;i<N;i++)
save[i]=0.0;
efft(save,N,0);//not add window
}
efft(x,N,windowflag);
y[0]=(float)(2*fabs(x[0])/N); //除N是为了使幅值与时域相等 smj 下同
y[N/2]=(float)(2*fabs(x[N/2])/N);
if(windowflag==2 && bRectify==0)// WINDOWFLAG=2 为 hanning window，note by smj
{
kt=0.0;
for(i=0;i<N;i++)
kt+=(float)(0.5*(1.0+cos(2*PI*i/N)));// 1/kt为函数窗的能量放大因子
for(i=1;i<N/2;i++)
y[i]=(float)(2*sqrt(x[i]*x[i]+x[N-i]*x[N-i])/kt/N);
}
else
{
for( i=1;i<N/2;i++)
y[i]=(float)(2*sqrt(x[i]*x[i]+x[N-i]*x[N-i])/N);//please note x/N
}
float maxoftotal=0.1; // Modified by XFY, 12/12/1999
short lastposiofmax=0; // Modified by XFY, 12/12/1999
for(i=0;i<spectrumpeaknumber;i++)
{
max =(float)( y[0]+1.0e-6); // Modified by XFY, 12/12/1999
for(j=0;j<=N/2;j++) // Modified by XFY, 12/12/1999
{
if(max<y[j])
{ max=y[j]; k=j;}
}
if(k==0) k=1;
if(i==0) // Modified by XFY, 12/12/1999
maxoftotal = (float)(y[k]+1.0e-6);
//if(bRectify) // Original
if(bRectify && (k-lastposiofmax) > 2 && y[k]>maxoftotal*0.1+1.0e-6) //Modified by XFY 12/12/1999
{
switch(windowflag)
{ case 1:
if(k<3||k>N/2-4)
break;
if(y[k+1]>=y[k-1])
dk=y[k+1]/(y[k]+y[k+1]); //dk 为频率校正公式 note by smj
else dk=-y[k-1]/(y[k]+y[k-1]);//参阅论文《频谱分析的校正方法》
amplititude=(float)(PI*dk*y[k]/sin(PI*dk));
frequence=(float)(samplefrequence*(k+dk)/N);
if(fabs(x[k])<=1e-6)
phase=(float)(PI/2);
else phase=(float)(atan(x[N-k]/x[k]));
if( x[k]<0.0 )
phase=phase+(float)PI;
if(x[k]>0.0 && x[N-k]<0.0)
phase=phase+(float)(2*PI);
phase=phase-(float)(dk*PI);
// phase=phase+(float)(dk*PI);
y[N-11*i-1]=amplititude;//校正幅值；smj
y[N-11*i-2]=frequence;
y[N-11*i-3]=phase;
break;
case 2:
if(k<3||k>N/2-4)
break;
if(y[k+1]>=y[k-1])
dk=(2*y[k+1]-y[k])/(y[k+1]+y[k]);
else dk=(y[k]-2*y[k-1])/(y[k]+y[k-1]);
amplititude=(float)(PI*dk*y[k]*2*(1-dk*dk)/sin(PI*dk));
frequence=samplefrequence*(k+dk)/N; //here I recommend dk, fre ammend is canceled.
//frequence=samplefrequence*k/N;
if(fabs(x[k])<1e-6)
phase=(float)PI/2;
else phase=(float)atan(x[N-k]/x[k]);
if( x[k]<-0.1 )
phase=phase+(float)PI;
if(x[k]>0.0 && x[N-k]<0.0)//smj recorrect 12/25/2000
// if(phase<-1.0e-005)
phase=phase+(float)(2*PI);
phase=phase-(float)(dk*PI);
// phase=phase+(float)(dk*PI);
y[N-11*i-1]=amplititude;
y[N-11*i-2]=frequence;
y[N-11*i-3]=phase;
break;
case 3:
case 4:
if(k<0||k>N/2-4)
break;
if(fabs(x[k])<1e-6)
ophase1=(float)PI/2;
else ophase1=(float)atan(x[N-k]/x[k]);
if(x[k]<0.0)
ophase1=ophase1+(float)PI;
if( x[k]>0.0 && x[N-k]<0.0 )
ophase1=ophase1+(float)(2*PI);
if(fabs(save[k])<1e-6)
ophase2=(float)(PI/2);
else ophase2=(float)atan(save[N-k]/save[k]);
if(save[k]<0.0)
ophase2=ophase2+(float)PI;
if( save[k]>0.0 && save[N-k]<0.0 )
ophase2=ophase2+(float)(2*PI);
dk=(float)(2*(ophase1-ophase2)/PI);
phase=ophase1-(float)(PI*dk*(N-1)/N);
if(windowflag==3)
amplititude=(float)(y[k]*PI/(dk*sin(PI*dk)*(0.54/(dk*dk)+0.46/(1-dk*dk))));
else
amplititude=(float)(y[k]*PI/(dk*sin(PI*dk)*(0.42/(dk*dk)+0.5/(1-dk*dk)-0.08/(4-dk*dk))));
frequence=samplefrequence*(k+dk)/N;
y[N-11*i-1]=amplititude;
y[N-11*i-2]=frequence;
y[N-11*i-3]=phase;
break;
}
}
else
{
y[N-11*i-1]=y[k];
y[N-11*i-2]=samplefrequence*k/N;
if(fabs(x[k])<1e-6)
phase=(float)PI/2;
else
phase=(float)atan(x[N-k]/x[k]);
if( x[k]<0.0 )
phase=phase+(float)PI;
if(x[k]>0.0 && x[N-k]<0.0)
phase=phase+(float)(2*PI);
y[N-11*i-3]=phase;
}
y[N-11*i-4]=(float)k;
if(k>=3 && k<=N/2-4)
{ for(j=k-2;j<=k+2;j++) //Modified by XFY, 16/12/1999
{ l=j-k+3;
y[N-11*i-5-l]=y[j];
y[j]=0.0;
}
}
else
{ y[N-11*i-5]=y[k];
if(k<3 && k!=0) //Modified by XFY, 12/12/1999
{ for(j=0;j<=k+2;j++)
y[j]=0.0;
}
else
{ for(j=k-2;j<=N/2;j++)
y[j]=0.0;
}
}
lastposiofmax = k; // Modified by XFY, 12/12/1999
}
for(i=0;i<spectrumpeaknumber;i++)
{ k=(int)y[N-11*i-4];
if(k>=3 && k<=N/2-4)
{ for(l=-2;l<=2;l++) //此处应为-2~+2 ？？？？？？？smj修改 4/23/2000
y[k+l]=y[N-i*11-(l+8)];
y[k]=y[N-11*i-1]; //amplititude has been amended
}
else
{
y[k]=y[N-11*i-5];
//y[N-11*i-5]=y[N-11*i-5]*kt;
}
}
for(i=0;i<N;i++)
x[i]=y[i];
//获取校正峰值，modified by XFY, 19/12/1999
for(i=0;i<spectrumpeaknumber;i++)
{
x[int(y[N-11*i-4])] = y[N-11*i-1];
// if(g_bSimulate && i==2 && windowflag>1)
// x[int(y[N-11*i-4])] = y[N-11*i-1]*2;
}
GlobalUnlock(hFloaty);
GlobalFree(hFloaty);
if(windowflag>2)
{ GlobalUnlock(hFloatsave);
GlobalFree(hFloatsave);
}
}
//////////趋势图用数据/////////////////////////////////////
// 趋势频谱提取
// x--输入实序列；
// n--序列长度；
// samplefrequence--采样频率；
// rpm--转速。
extern "C" __declspec(dllexport)void amendfft1( float *x, int n, float samplefrequence,float rpm )
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,j,k,mm,position;
float max,phase[513];
float dk,aver=0.0;
//原始数据去均值2001/04/09
for(i=0;i<n;i++)
aver+=x[i];
aver=aver/n;
for(i=0;i<n;i++)
x[i]-=aver;
// m = int (log(float(n))/log(2.0));
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?