📄 dph.c
字号:
#include <windows.h>
#include <math.h>
#define PI 3.141592654
struct sucoeff{
double ug0;
double ut0;
} ;
struct sucoeff usinbal ;//定义结构体usinbal为由单面平衡所得的系统的影响系数(复数），用模ug0 与 幅角ut0(rad)表示
struct sucoeff gm; //定义结构体gm为平衡面上的原始不平衡量，幅角单位为度。
struct sucoeff gi; //定义结构体gi为平衡面上由试重（不可拆）产生的新增不平衡量，幅角单位为度。
struct sucoeff phres; //定义结构体phres为平衡面上应加配重的大小与相位，相位单位为度。
struct strsb {
double g0;
double t0;
struct sucoeff usinbal;
} ;
struct strsb sinbal; //定义结构体sinbal为由单面平衡所得的不平衡量的大小（g）和相位(角度）以及影响系数。
struct ducoeff{
double du11;
double dthe11;
double du12;
double dthe12;
double du21;
double dthe21;
double du22;
double dthe22;
};
struct ducoeff udoubal; //定义结构体udoubal为由双面平衡所得的系统的影响系数(4个复数），
// 用模du11 与 幅角dth11(rad)等表示;
struct strdb{
double g10;
double g20;
double t10;
double t20;
struct ducoeff udoubal;
};
struct strdb doubal; //定义结构体doubal为由双面平衡所得的不平衡量的大小(g)和相位（角度）
struct strcalf{
double a;
double b;
} ;
struct strcalf aphase; //DFT算法的结果（复数），代表振动值，是全局变量。
//////////////////////////////////////////////////////////////////////////////////////
struct strcalf FAR PASCAL calf(double ab[]);
struct strsb FAR PASCAL sb(double sa0[],double sa1[],double g,double bata);
struct strsb FAR PASCAL sonly(double sa[],struct sucoeff usinbal);
struct strdb FAR PASCAL db(double da10[],double da20[],double da11[],double da12[],double da21[],double da22[],
double g1,double bata1,double g2,double bata2);
struct strdb FAR PASCAL donly(double da10[],double da20[],struct ducoeff udoubal);
void FAR PASCAL fft(double pr[],double pi[],int n,int k,double fr[],double fi[]);
struct sucoeff FAR PASCAL ph(struct sucoeff gm,struct sucoeff gi);
//////////////////////////////////////////////////////////////////////////////////////
int FAR PASCAL LibMain(HANDLE hInst,WORD wData,WORD wHeap,LPSTR lpCmd)
{
if(wHeap) UnlockData(0);
return 1;
}
int FAR PASCAL WEP(int nSysExt)
{
return 1;
}
////////////////////////////////////////////////////////////
//单面平衡算法：采用加一次配重方法，求系统的原始不平衡量，
// 同时获取所测系统的影响系数。
struct strsb FAR PASCAL sb(double sa0[32],double sa1[32],double g,double bata)
{
double a0,afa0,a1,afa1;
double fz,fm,u,b,c;
double g0,t0;
calf(sa0);
a0=aphase.a;
afa0=aphase.b;
calf(sa1);
a1=aphase.a;
afa1=aphase.b;
fz=a1*sin(afa1)-a0*sin(afa0);
fm=a1*cos(afa1)-a0*cos(afa0);
u=sqrt(fm*fm+fz*fz);
u=u/g;
b=atan(fz/fm);
if (fm<0.0)
b=b+PI;
if (b<0.0)
b=b+2.0*PI;
c=b-bata*PI/180.0;
if (c<0.0)
c=c+2.0*PI;
if (c>2.0*PI)
c=c-2.0*PI;
sinbal.usinbal.ug0=u;
sinbal.usinbal.ut0=c;
g0=a0/u;
t0=afa0-c;
t0=t0*180.0/PI;
if (t0<0.0)
t0=t0+360.0;
if (t0>360.0)
t0=t0-360.0;
sinbal.g0=g0;
sinbal.t0=t0;
return sinbal;
}
//单面平衡算法：用于实时监测系统的残余不平衡量。
//说明：已知系统的影响系数 usinbal，振动幅值 sa[32] 已测得，该函数可计算出
// 系统的残余不平衡量，存放于结构体 sinbal 中。
struct strsb FAR PASCAL sonly(double sa[32],struct sucoeff usinbal)
{
double a0,afa0,bug0,but0;
double g0,t0;
calf(sa);
a0=aphase.a;
afa0=aphase.b;
bug0=usinbal.ug0;
but0=usinbal.ut0;
g0=a0/bug0;
t0=afa0-but0;
t0=t0*180.0/PI;
if (t0<0.0)
t0=t0+360.0;
if (t0>360.0)
t0=t0-360.0;
sinbal.g0=g0;
sinbal.t0=t0;
return sinbal;
}
//////////////////////////////////////////////////////////////////////
//双面平衡算法：采用加两次配重法，求系统的原始不平衡量，
// 同时获取所测系统的影响系数。
struct strdb FAR PASCAL db(double da10[32],double da20[32],double da11[32],double da12[32],double da21[32],
double da22[32],double g1,double bata1,double g2,double bata2)
{
double a10,t10,a20,t20,a11,t11,a12,t12,a21,t21,a22,t22;
double u11,the11,u12,the12,u21,the21,u22,the22;
double ca11,ca12,ca13,ca14;
double ca21,ca22,ca23,ca24;
double ca31,ca32,ca33,ca34;
double ca41,ca42,ca43,ca44;
double c11,c12,c13;
double c21,c22,c23;
double c31,c32,c33;
double d1,d2,d3,b1,b2,b3,b4,e11,e12,e21,e22,f1,f2;
double x1,x2,x3,x4;
double ur11,ui11,ur12,ui12,ur21,ui21,ur22,ui22;
calf(da10);
a10=aphase.a;
t10=aphase.b;
calf(da20);
a20=aphase.a;
t20=aphase.b;
calf(da11);
a11=aphase.a;
t11=aphase.b;
calf(da12);
a12=aphase.a;
t12=aphase.b;
calf(da21);
a21=aphase.a;
t21=aphase.b;
calf(da22);
a22=aphase.a;
t22=aphase.b;
u11=sqrt(pow((a11*cos(t11)-a10*cos(t10)),2)+pow((a11*sin(t11)-a10*sin(t10)),2));
u11=u11/g1;
u12=sqrt(pow((a12*cos(t12)-a11*cos(t11)),2)+pow((a12*sin(t12)-a11*sin(t11)),2));
u12=u12/g2;
u21=sqrt(pow((a21*cos(t21)-a20*cos(t20)),2)+pow((a21*sin(t21)-a20*sin(t20)),2));
u21=u21/g1;
u22=sqrt(pow((a22*cos(t22)-a21*cos(t21)),2)+pow((a22*sin(t22)-a21*sin(t21)),2));
u22=u22/g2;
doubal.udoubal.du11=u11;
doubal.udoubal.du12=u12;
doubal.udoubal.du21=u21;
doubal.udoubal.du22=u22;
the11=atan((a11*sin(t11)-a10*sin(t10))/(a11*cos(t11)-a10*cos(t10)));
the12=atan((a12*sin(t12)-a11*sin(t11))/(a12*cos(t12)-a11*cos(t11)));
the21=atan((a21*sin(t21)-a20*sin(t20))/(a21*cos(t21)-a20*cos(t20)));
the22=atan((a22*sin(t22)-a21*sin(t21))/(a22*cos(t22)-a21*cos(t21)));
if ((a11*cos(t11)-a10*cos(t10))<0.0)
the11=the11+PI;
if ((a12*cos(t12)-a11*cos(t11))<0.0)
the12=the12+PI;
if ((a21*cos(t21)-a20*cos(t20))<0.0)
the21=the21+PI;
if ((a22*cos(t22)-a21*cos(t21))<0.0)
the22=the22+PI;
if (the11<0.0)
the11=the11+2.0*PI;
if (the12<0.0)
the12=the12+2.0*PI;
if (the21<0.0)
the21=the21+2.0*PI;
if (the22<0.0)
the22=the22+2.0*PI;
the11=the11-bata1*PI/180.0;
the21=the21-bata1*PI/180.0;
the12=the12-bata2*PI/180.0;
the22=the22-bata2*PI/180.0;
doubal.udoubal.dthe11=the11;
doubal.udoubal.dthe21=the21;
doubal.udoubal.dthe12=the12;
doubal.udoubal.dthe22=the22;
ur11=u11*cos(the11);
ui11=u11*sin(the11);
ur12=u12*cos(the12);
ui12=u12*sin(the12);
ur21=u21*cos(the21);
ui21=u21*sin(the21);
ur22=u22*cos(the22);
ui22=u22*sin(the22);
ca11=ur11;
ca12=-ui11;
ca13=ur12;
ca14=-ui12;
ca21=ui11;
ca22=ur11;
ca23=ui12;
ca24=ur12;
ca31=ur21;
ca32=-ui21;
ca33=ur22;
ca34=-ui22;
ca41=ui21;
ca42=ur21;
ca43=ui22;
ca44=ur22;
b1=a10*cos(t10);
b2=a10*sin(t10);
b3=a20*cos(t20);
b4=a20*sin(t20);
c11=ca21*ca12-ca11*ca22;
c12=ca21*ca13-ca11*ca23;
c13=ca21*ca14-ca11*ca24;
c21=ca31*ca12-ca11*ca32;
c22=ca31*ca13-ca11*ca33;
c23=ca31*ca14-ca11*ca34;
c31=ca41*ca12-ca11*ca42;
c32=ca41*ca13-ca11*ca43;
c33=ca41*ca14-ca11*ca44;
d1=ca21*b1-ca11*b2;
d2=ca31*b1-ca11*b3;
d3=ca41*b1-ca11*b4;
e11=c21*c12-c11*c22;
e12=c21*c13-c11*c23;
e21=c31*c12-c11*c32;
e22=c31*c13-c11*c33;
f1=c21*d1-c11*d2;
f2=c31*d1-c11*d3;
x4=(e21*f1-e11*f2)/(e21*e12-e11*e22);
x3=(f1-e12*x4)/e11;
x2=(d1-c12*x3-c13*x4)/c11;
x1=(b1-ca12*x2-ca13*x3-ca14*x4)/ca11;
a10=sqrt(x1*x1+x2*x2);
a20=sqrt(x3*x3+x4*x4);
t10=atan(x2/x1);
t20=atan(x4/x3);
if (x1<0.0)
t10=t10+PI;
if (x3<0.0)
t20=t20+PI;
t10=t10*180.0/PI;
t20=t20*180.0/PI;
if (t10<0.0)
t10=t10+360.0;
if (t10>360.0)
t10=t10-360.0;
if (t20<0.0)
t20=t20+360.0;
if (t20>360.0)
t20=t20-360.0;
doubal.g10=a10;
doubal.g20=a20;
doubal.t10=t10;
doubal.t20=t20;
return doubal;
}
//双面平衡算法：用于实时监测系统的残余不平衡量。
//说明：已知系统的影响系数 udoubal，第一支承处和第二支承处的振动幅值 da10[32] ,da20[32]已测得，
// 该函数可计算出系统第一平衡面和第二平衡面上的残余不平衡量，结果存放于结构体 doubal 中。
struct strdb FAR PASCAL donly(double da10[32],double da20[32],struct ducoeff udoubal)
{
double a10,t10,a20,t20;
double u11,the11,u12,the12,u21,the21,u22,the22;
double ca11,ca12,ca13,ca14;
double ca21,ca22,ca23,ca24;
double ca31,ca32,ca33,ca34;
double ca41,ca42,ca43,ca44;
double c11,c12,c13;
double c21,c22,c23;
double c31,c32,c33;
double d1,d2,d3,b1,b2,b3,b4,e11,e12,e21,e22,f1,f2;
double x1,x2,x3,x4;
double ur11,ui11,ur12,ui12,ur21,ui21,ur22,ui22;
u11=udoubal.du11;
u12=udoubal.du12;
u21=udoubal.du21;
u22=udoubal.du22;
the11=udoubal.dthe11;
the21=udoubal.dthe21;
the12=udoubal.dthe12;
the22=udoubal.dthe22;
calf(da10);
a10=aphase.a;
t10=aphase.b;
calf(da20);
a20=aphase.a;
t20=aphase.b;
ur11=u11*cos(the11);
ui11=u11*sin(the11);
ur12=u12*cos(the12);
ui12=u12*sin(the12);
ur21=u21*cos(the21);
ui21=u21*sin(the21);
ur22=u22*cos(the22);
ui22=u22*sin(the22);
ca11=ur11;
ca12=-ui11;
ca13=ur12;
ca14=-ui12;
ca21=ui11;
ca22=ur11;
ca23=ui12;
ca24=ur12;
ca31=ur21;
ca32=-ui21;
ca33=ur22;
ca34=-ui22;
ca41=ui21;
ca42=ur21;
ca43=ui22;
ca44=ur22;
b1=a10*cos(t10);
b2=a10*sin(t10);
b3=a20*cos(t20);
b4=a20*sin(t20);
c11=ca21*ca12-ca11*ca22;
c12=ca21*ca13-ca11*ca23;
c13=ca21*ca14-ca11*ca24;
c21=ca31*ca12-ca11*ca32;
c22=ca31*ca13-ca11*ca33;
c23=ca31*ca14-ca11*ca34;
c31=ca41*ca12-ca11*ca42;
c32=ca41*ca13-ca11*ca43;
c33=ca41*ca14-ca11*ca44;
d1=ca21*b1-ca11*b2;
d2=ca31*b1-ca11*b3;
d3=ca41*b1-ca11*b4;
e11=c21*c12-c11*c22;
e12=c21*c13-c11*c23;
e21=c31*c12-c11*c32;
e22=c31*c13-c11*c33;
f1=c21*d1-c11*d2;
f2=c31*d1-c11*d3;
x4=(e21*f1-e11*f2)/(e21*e12-e11*e22);
x3=(f1-e12*x4)/e11;
x2=(d1-c12*x3-c13*x4)/c11;
x1=(b1-ca12*x2-ca13*x3-ca14*x4)/ca11;
a10=sqrt(x1*x1+x2*x2);
a20=sqrt(x3*x3+x4*x4);
t10=atan(x2/x1);
t20=atan(x4/x3);
if (x1<0.0)
t10=t10+PI;
if (x3<0.0)
t20=t20+PI;
t10=t10*180.0/PI;
t20=t20*180.0/PI;
if (t10<0.0)
t10=t10+360.0;
if (t10>360.0)
t10=t10-360.0;
if (t20<0.0)
t20=t20+360.0;
if (t20>360.0)
t20=t20-360.0;
doubal.g10=a10;
doubal.g20=a20;
doubal.t10=t10;
doubal.t20=t20;
return doubal;
}
//////////////////////////////////////////////////////////
//DFT算法：
struct strcalf FAR PASCAL calf(double ab[32])
{
double a1,b1,a,b;
int i;
a1=0.0;
b1=0.0;
for (i=0;i<32;i=i+1) {
a1=a1+ab[i]*cos(i*PI/16.0);
b1=b1-ab[i]*sin(i*PI/16.0);
};
a=sqrt(a1*a1+b1*b1)/16.0;
b=atan(b1/a1);
if (a1<0.0)
b=b+PI;
aphase.a=a;
aphase.b=b;
return aphase;
}
//////////////////////////////////////////////////////////////
//FFT算法
//注：数组pr[]，pi[]（令pi[]=0）为时域采样信号值，n为采样点数，K=lg（n）；
// FFT结果有两种表现形式：数组pr[]，pi[]以模与幅角(单位度)的形式存放FFT结果；
// 数组fr[]，fi[]以实部与虚部的形式存放FFT结果；
void FAR PASCAL fft(double pr[],double pi[],int n,int k,double fr[],double fi[])
{ 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);
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;
}
}
for (i=0; i<=n-1; i++)
{ pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if (fabs(fr[i])<0.000001*fabs(fi[i]))
{ if ((fi[i]*fr[i])>0) pi[i]=90.0;
else pi[i]=-90.0;
}
else
pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
}
return;
}
//////////////////////////////////////////////////////////////
//配重算法ph：考虑到试重不可拆，ph函数解决平衡面上实际应加配重的大小与相位（单位：度）；
//注： 定义结构体gm为平衡面上的原始不平衡量，幅角单位为度。
// 定义结构体gi为平衡面上由试重（不可拆）产生的新增不平衡量，幅角单位为度。
// 定义结构体phres为平衡面上应加配重的大小与相位，相位单位为度。
struct sucoeff FAR PASCAL ph(struct sucoeff gm,struct sucoeff gi)
{
struct sucoeff phres;
double gx,gy,tp;
double g,t,g1,bata1;
g=gm.ug0;
t=gm.ut0*PI/180.0;
g1=gi.ug0;
bata1=gi.ut0*PI/180.0;
gx=g*cos(t)+g1*cos(bata1);
gy=g*sin(t)+g1*sin(bata1);
phres.ug0=sqrt(gx*gx+gy*gy);
if(fabs(gx)<0.000001*fabs(gy))
{
if((gx*gy)>0)
tp=90.0;
else
tp=-90.0;
}else
tp=atan(gy/gx)*360.0/6.283185306;
if(gx<0)
tp=tp+180.0;
if(tp<180.0)
{
tp=tp+180.0;
}else
tp=tp-180.0;
phres.ut0=tp;
return phres;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -