signalcalc.cpp
来自「常用算法与数据结构原代码」· C++ 代码 · 共 1,624 行 · 第 1/3 页
CPP
1,624 行
efft(x,n,2);//window 2
//if(x[0]<0) phase[0]=180.0;
//else phase[0]=0.0;
if(x[n/2]<0) phase[n/2]=180.0;
else phase[n/2]=0.0;
for(i=1;i<n/2;i++)
{ if(fabs(x[i])<1.0e-6) phase[i] = (float)PI/2;
else phase[i] = (float)atan(x[n-i]/x[i]);
phase[i]=(float)(phase[i]*180./PI);
if(x[i]<0) phase[i]=180+phase[i];
else if((x[i]>0)&&(x[n-i]<0)) phase[i]=360+phase[i];
}
x[0] = 0.0; // for mm=1; while mm=1, mm / 2 = 0;
phase[0] = 0.0;
x[n/2]=(float)(2*fabs(x[n/2])/n);
for( i=1;i<n/2;i++)
x[i]=(float)(2*sqrt(x[i]*x[i]+x[n-i]*x[n-i])/n);//please note x/n
mm=(int)(rpm/60./(samplefrequence/n));//basic fre
if(mm <= 5) //smj modified,2001/2/27
{
for(i=1;i<=5;i++)
{
switch(i)
{
case 1: position=(int)mm/2; break;
case 2: position=mm; break;
case 3: position=2*mm; break;
case 4: position=3*mm; break;
case 5: position=4*mm; break;
}
x[n-i] = x[position];
if(phase[position]>0) phase[position] = (int)phase[position] % 360;
if(phase[position]<0) phase[position] = (int)phase[position] % 360 + 360;
x[n-i-8] = phase[position];
}
return;
}
int lastposi = 0;
int searchrange=2;//smj 2001/2/27
/*
if (mm>25) searchrange=6;
else if (mm>21) searchrange=5;
else if (mm>17) searchrange=4;
else if (mm>13) searchrange=3;
else if (mm>9) searchrange=2;
else searchrange=1;
*/
for(i=1;i<=5;i++)
{
switch(i)
{ case 1: position=(int)mm/2; break;
case 2: position=mm; break;
case 3: position=2*mm; break;
case 4: position=3*mm; break;
case 5: position=4*mm; break;
}
max=-900.0;
for(j=position-searchrange; j<=position+searchrange; j++)//search range,modified by smj 2001/2/27
if(j>=0 && max<x[j])
{ max=x[j];
k=j;
}
if(k>=3 && k<=n/2-4 && x[k] != 0.0 && k != lastposi)
{
if(x[k+1]>=x[k-1])
dk=(2*x[k+1]-x[k])/(x[k+1]+x[k]);
else
dk=(x[k]-2*x[k-1])/(x[k]+x[k-1]);
x[k]=(float)(PI*dk*x[k]*2*(1-dk*dk)/sin(PI*dk));
phase[k]=(float)(phase[k]-dk*180);
if(phase[k]>360) phase[k]=phase[k]-360;
if(phase[k]<0) phase[k]=phase[k]+360;
}
lastposi = k;
x[n-i] = x[k];
if(phase[k]>0) phase[k] = (int)phase[k] % 360;
if(phase[k]<0) phase[k] = (int)phase[k] % 360 + 360;
x[n-i-8] = phase[k];
}
}
//////////趋势图用数据/////////////////////////////////////
extern "C" __declspec(dllexport)void amendfft2( float *x, int n, float samplefrequence,float *FreqTable, int freqnum )
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,j,k,position;
float max,phase[513];
float 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));
efft(x,n,2);//window 2
//if(x[0]<0) phase[0]=180.0;
//else phase[0]=0.0;
if(x[n/2]<0) phase[n/2]=180.0;
else phase[n/2]=0.0;
for(i=1;i<n/2;i++)
{ if(fabs(x[i])<1.0e-6) phase[i] = (float)PI/2;
else phase[i] = (float)atan(x[n-i]/x[i]);
phase[i]=(float)(phase[i]*180./PI);
if(x[i]<0) phase[i]=180+phase[i];
else if((x[i]>0)&&(x[n-i]<0)) phase[i]=360+phase[i];
}
x[0] = 0.0; // for mm=1; while mm=1, mm / 2 = 0;
phase[0] = 0.0;
x[n/2]=(float)(2*fabs(x[n/2])/n);
for( i=1;i<n/2;i++)
x[i]=(float)(2*sqrt(x[i]*x[i]+x[n-i]*x[n-i])/n);//please note x/n
// mm=(int)(rpm/60./(samplefrequence/n));//basic fre
for(i=1; i<=(freqnum>8 ? 8: freqnum);i++)
{
position = (int)FreqTable[i-1]/(samplefrequence/n);
max=-900.0;
for(j=position-2; j<=position+2; j++)//search range,modified by smj 2001/2/27
if(j>=0 && max<x[j])
{ max=x[j];
k=j;
}
x[n-i] = max;
if(phase[k]>0) phase[k] = (int)phase[k] % 360;
if(phase[k]<0) phase[k] = (int)phase[k] % 360 + 360;
x[n-i-8] = phase[k];
}
/*
if(mm <= 5) //smj modified,2001/2/27
{
for(i=1;i<=5;i++)
{
switch(i)
{
case 1: position=(int)mm/2; break;
case 2: position=mm; break;
case 3: position=2*mm; break;
case 4: position=3*mm; break;
case 5: position=4*mm; break;
}
x[n-i] = x[position];
if(phase[position]>0) phase[position] = (int)phase[position] % 360;
if(phase[position]<0) phase[position] = (int)phase[position] % 360 + 360;
x[n-i-8] = phase[position];
}
return;
}
int lastposi = 0;
int searchrange;//smj 2001/2/27
if (mm>25) searchrange=6;
else if (mm>21) searchrange=5;
else if (mm>17) searchrange=4;
else if (mm>13) searchrange=3;
else if (mm>9) searchrange=2;
else searchrange=1;
for(i=1;i<=5;i++)
{
switch(i)
{ case 1: position=(int)mm/2; break;
case 2: position=mm; break;
case 3: position=2*mm; break;
case 4: position=3*mm; break;
case 5: position=4*mm; break;
}
max=-900.0;
for(j=position-searchrange; j<=position+searchrange; j++)//search range,modified by smj 2001/2/27
if(j>=0 && max<x[j])
{ max=x[j];
k=j;
}
if(k>=3 && k<=n/2-4 && x[k] != 0.0 && k != lastposi)
{
if(x[k+1]>=x[k-1])
dk=(2*x[k+1]-x[k])/(x[k+1]+x[k]);
else
dk=(x[k]-2*x[k-1])/(x[k]+x[k-1]);
x[k]=(float)(PI*dk*x[k]*2*(1-dk*dk)/sin(PI*dk));
phase[k]=(float)(phase[k]-dk*180);
if(phase[k]>360) phase[k]=phase[k]-360;
if(phase[k]<0) phase[k]=phase[k]+360;
}
lastposi = k;
x[n-i] = x[k];
if(phase[k]>0) phase[k] = (int)phase[k] % 360;
if(phase[k]<0) phase[k] = (int)phase[k] % 360 + 360;
x[n-i-8] = phase[k];
}*/
}
////////////////////倒谱计算///////////////////////////////////////////
// 幅值倒谱计算
// x--输入实序列，输出的前N/2为倒谱；
// N--序列长度；
extern "C" __declspec(dllexport)void rceps(float *xr,int N)// 参考MATLAB得到的幅值倒谱 smj
{ //MATLAB中的算法为xhat = real(ifft(log(abs(fft(x)))))
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,m0;
float *xi;
HANDLE hFloatxi;
hFloatxi=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
xi=(float *)GlobalLock(hFloatxi);
//原始数据去均值2001/03/06
float aver=0.0f;
for(i=0;i<N;i++) aver+=xr[i];
aver=aver/N;
for(i=0;i<N;i++) xr[i]-=aver;
for(i=0;i<N;i++)
xi[i]=0.0f;
m0=CalculatePower(N);
fft( xr,xi,m0,0 );
for(i=0;i<N;i++)
{
xr[i]=(float)sqrt(xr[i]*xr[i]+xi[i]*xi[i]);
if(xr[i]==0)xr[i]=1.0e-6;
xr[i]=(float)log(double(xr[i]));
xi[i]=0.0f;
}
fft(xr,xi,m0,1);
GlobalUnlock(hFloatxi);
GlobalFree(hFloatxi);
}
extern "C" __declspec(dllexport) void waterfallfft( float *x,int n,float samplefrequence,int windowflag)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,k;
float max,amplitude,phase,frequence;
float dk;
max=-10.0;
float* y = NULL;
y = (float *)::VirtualAlloc(NULL, sizeof(float) * (n/2+1),
MEM_RESERVE | MEM_COMMIT, PAGE_READWRITE);
::ZeroMemory(y, sizeof(float) * n/2);
efft(x,n,windowflag);
y[0]=(float)(2*fabs(x[0])/n);
y[n/2]=(float)(2*fabs(x[n/2])/n);
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
for(i=29;i<35;i++)
{ if(max<y[i])
{ max=y[i];
k=i;
}
}
switch(windowflag)
{ case 0:
if(y[k+1]>=y[k-1])
dk=y[k+1]/(y[k]+y[k+1]);
else dk=-y[k-1]/(y[k]+y[k-1]);
amplitude=(float)(PI*dk*y[k]/sin(PI*dk));
frequence=samplefrequence*(k+dk)/n;
phase=(float)atan(x[n-k]/x[k]);
if( x[n-k]<0.0 && x[k]<0.0)
phase=phase-(float)PI;
if( x[k]<0.0 && x[n-k]>0.0 )
phase=phase+(float)PI;
phase=phase-(float)(dk*PI);
break;
case 2:
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]);
amplitude=(float)(PI*dk*y[k]*2*(1-dk*dk)/sin(PI*dk));
frequence=samplefrequence*(k+dk)/n;
phase=(float)atan(x[n-k]/x[k]);
if( x[n-k]<0.0 && x[k]<0.0)
phase=phase-(float)PI;
if( x[k]<0.0 && x[n-k]>0.0 )
phase=phase+(float)PI;
phase=phase-(float)(dk*PI); //phase in -180-180
break;
}
if(phase>2*PI)
phase=phase-(float)(2*PI);
if(phase<0)
phase=phase+(float)(2*PI);
phase=(float)(phase*180.0/PI);
x[n-1]=frequence;
x[n-2]=amplitude;
x[n-3]=phase;
for(i=0;i<n/2;i++)
x[i]=y[i];
::VirtualFree(y, sizeof(float) * (n/2+1), MEM_DECOMMIT);
::VirtualFree(y, 0, MEM_RELEASE);
}
//////////////////////希尔伯特变换及包络////////////////////////////////////////
extern "C" __declspec(dllexport)void Hilbert(float *x,float *xi,int N)
{ //N为信号序列长度，
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int B,i;
float *y;
y=(float *)malloc(N*sizeof(float));
// xi=(float *)malloc(N*sizeof(float));
for(i=0;i<N;i++)
{
y[i]=x[i];
xi[i]=0.0f;
}
B=CalculatePower(N);
fft(y,xi,B,0);
for(i=1;i<N/2;i++)
{
y[i]=2*y[i];
xi[i]=2*xi[i];
}
for(i=N/2;i<N;i++)
{
y[i]=0.0f;
xi[i]=0.0f;
}
fft(y,xi,B,1);
free(y);
}
/*
extern "C" __declspec(dllexport)void Hilbert(float *x,float *xi,int N)
{ //N为信号序列长度，
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int B,i;
float *y;
y=(float *)malloc(N*sizeof(float));
// xi=(float *)malloc(N*sizeof(float));
for(i=0;i<N;i++)
{
y[i]=x[i];
xi[i]=0.0f;
}
B=CalculatePower(N);
fft(y,xi,B,0);
for(i=1;i<N/2;i++)
{
y[i]=2*y[i];
xi[i]=2*xi[i];
}
for(i=N/2;i<N;i++)
{
y[i]=0.0f;
xi[i]=0.0f;
}
fft(y,xi,B,1);
free(y);
}*/
// 计算信号的希尔伯特包络
// x--输入实序列，输出为包络；
// N--序列长度。
extern "C" __declspec(dllexport)void Envelope(float *x,int N)
{ //N为信号序列长度，
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i;
float *xi;
xi=(float *)malloc(N*sizeof(float));
for(i=0;i<N;i++)
xi[i]=0.0f;
Hilbert(x,xi,N);
for(i=0;i<N;i++)
x[i]=(float)sqrt(xi[i]*xi[i]+x[i]*x[i]);
}
//////////////////////////相关函数/////////////////////////////////////////
// x--输入实序列，输出为相关函数；
// y--输入实序列；
// N--序列长度；
extern "C" __declspec(dllexport) void correlation(float *x,float *y ,int N)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,j;
float sumx,sumy,averx,avery;
float *z;
HANDLE hFloatz;
hFloatz=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,(2*N-1)*sizeof(float));
z=(float *)GlobalLock(hFloatz);
sumx=sumy=0.0;
for(i=0;i<N;i++)
{
sumx+=x[i];
sumy+=y[i];
}
averx=sumx/N;
avery=sumy/N;
for(i=0;i<N;i++)
{
x[i]-=averx;
y[i]-=avery;
}
for(j=-N+1;j<N;j++)
{
z[j+N-1]=0.0;
for(i=0;i<N;i++)
{
if(((i+j)>=0)&&((i+j)<N))
z[j+N-1]+=x[i]*y[i+j];
}
}
for(i=0;i<2*N-1;i++)
x[i]=z[i]/N;
GlobalUnlock(hFloatz);
GlobalFree(hFloatz);
}
///////////////////////功率谱计算////////////////////////////////////
// 互谱估计
// x-输入序列，输出的前nfft/2+1个值为互功率谱值；
// y-输入序列；
// N-输入序列的长度，为2的幂次方；
// M-分段的长度，为2的幂次方；(M<=N)
// lag-估计功率谱所用相关函数的点数，其值不大于M/2+1；
// nfft-估计功率谱所用FFT长度，不小于2lag-1；
// windowflag- 窗函数类型；
extern "C" __declspec(dllexport)void csd(float *x,float *y,int N,int M,int lag,int nfft,int windowflag)
{ //N为信号序列长度，
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,j,k,m,nSect;
float sumx,sumy,averx,avery;
float *u,*v,*z,*zi;
HANDLE hFloatu,hFloatv,hFloatz,hFloatzi;
hFloatu=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,M*sizeof(float));
u=(float *)GlobalLock(hFloatu);
hFloatv=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,M*sizeof(float));
v=(float *)GlobalLock(hFloatv);
hFloatz=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
z=(float *)GlobalLock(hFloatz);
hFloatzi=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
zi=(float *)GlobalLock(hFloatzi);
sumx=sumy=0.0;
for(i=0;i<N;i++)
{
sumx+=x[i];
sumy+=y[i];
}
averx=sumx/N;
avery=sumy/N;
for(i=0;i<N;i++)
{
x[i]=x[i]-averx;
y[i]=y[i]-avery;
}
nSect=(N-M/2)/(M/2);
for(i=0;i<M/2+1;i++)
{
z[i]=0.0;
}
for(k=0;k<nSect;k++)
{
for(i=0;i<M;i++)
{
u[i]=x[k*M/2+i];
}
for(i=0;i<M/2;i++)
{
v[i]=y[k*M/2+i];
}
for(i=M/2;i<M;i++)
{
v[i]=0.0;
}
for(i=0;i<M/2+1;i++)
{
zi[i]=0.0;
for(j=0;j<M;j++)
{
if(((i+j)>=0)&&((i+j)<M))
zi[i]+=u[j+i]*v[j]; //zi存放每段的互相关值zi[2]=x2y0+x3y1+...+x(1+M/2)y(M/2-1);
}
}
for(j=0;j<M/2+1;j++)
z[j]+=zi[j]; //z将每段的第 j个互相关值累加；
}
for(i=0;i<M/2+1;i++)
{
z[i]=z[i]/N;
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?