signalcalc.cpp
来自「常用算法与数据结构原代码」· C++ 代码 · 共 1,624 行 · 第 1/3 页
CPP
1,624 行
}
windowed( z,lag,windowflag);
for(i=0;i<nfft;i++)
zi[i]=0.0;
for(i=lag;i<nfft;i++)
z[i]=0.0;
m=CalculatePower(nfft);
fft(z,zi,m,0);
for(i=0;i<nfft/2+1;i++)
x[i]=(float)(20.0*log10(sqrt(z[i]*z[i]+zi[i]*zi[i])));
GlobalUnlock(hFloatu);
GlobalUnlock(hFloatv);
GlobalUnlock(hFloatz);
GlobalUnlock(hFloatzi);
GlobalFree(hFloatu);
GlobalFree(hFloatv);
GlobalFree(hFloatz);
GlobalFree(hFloatzi);
}
///////////////////离散小波变换///////////////////////////////////////
int flag[max_nscales+1],sca[max_nscales];
/* 离散小波变换*///采用Mallat塔式算法 smj
/* a存放原信号与小波分解的各级逼近信号；
d存放小波分解的各级细节信号,长度为2N-N/pow(2,J)
wlen为小波长度
nscales为小波分解的级数
g为尺度系数
h为小波系数*/
void mallat(float *a,float *d,float *h,float *g,int wlen,int nscales)
{
int i,j,k,mid;
float p,q;
for(j=1;j<=nscales;j++)//nscales为小波展开的级数的大小
{
for(i=0;i<sca[j];i++)
{
p=0.0;
q=0.0;
for(k=0;k<wlen;k++)
{
mid=k+2*i;
if(mid>=sca[j-1]) mid=mid-sca[j-1];
p=p+h[k]*a[flag[j-1]+mid];
q=q+g[k]*a[flag[j-1]+mid];
}
a[flag[j]+i]=p;
d[flag[j]+i]=q;
}
}
// free(flag);
}
/* 离散小波反变换*///采用Mallat塔式算法 smj
/* a的输出为重建信号 */
void imallat(float *a,float *d,float *h,float *g,int wlen,int nscales)
{
int i,j,k,mid;
float p,q;
/* flag=(int*)malloc(nscales*sizeof(int));
flag[0]=0;
for(i=0;i<nscales;i++)
flag[i+1]=flag[i]+sca[i];*/
for(k=nscales;k>0;k--)//nscales为小波展开的级数的大小
{
for(i=0;i<sca[k];i++)
{
p=0.0;
q=0.0;
for(j=0;j<wlen/2;j++)
{
mid=i-j;
if(mid<0) mid=sca[k]+(i-j);
p+=h[2*j]*a[flag[k]+mid]+g[2*j]*d[flag[k]+mid];
q+=h[2*j+1]*a[flag[k]+mid]+g[2*j+1]*d[flag[k]+mid];
}
a[flag[k-1]+2*i]=p;
a[flag[k-1]+2*i+1]=q;
}
}
// free(flag);
}
// 离散小波变换
// a-存放原信号与小波分解的各级逼近信号；
// d-存放小波分解的各级细节信号,长度为2N-N/pow(2,J)
// N-原信号长度
// filter_kind-小波类型3~10分别为db3~db10小波
// nscales-小波分解的级数
// inv-0为正变换，1为逆变换
extern "C" __declspec(dllexport) void DWT(float *a,float *d,int N,int filter_kind,int nscales,int inv)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int i,j,wlen;
float h[max_wlen],g[max_wlen];
switch(filter_kind)
{
case 3:
wlen=6;//Daubiches3小波
h[0]=0.332670552950;
h[1]=0.806891509311;
h[2]=0.459877502118;
h[3]=-0.135011020010;
h[4]=-0.085441273882;
h[5]=0.035226291882;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 4:
wlen=8;//Daubiches4小波
h[0]=0.2303778133089;
h[1]=0.714846570553;
h[2]=0.630880767930;
h[3]=-0.027983769417;
h[4]=-0.187034811719;
h[5]=0.030841381836;
h[6]=0.032883011667;
h[7]=-0.010597401785;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 5:
wlen=10;//Daubiches5小波
h[0]=0.160102397974;
h[1]=0.603829269797;
h[2]=0.724308528438;
h[3]=0.138428145901;
h[4]=-0.242294887066;
h[5]=-0.032244869585;
h[6]=0.077571493840;
h[7]=-0.006241490213;
h[8]=-0.012580751999;
h[9]=0.003335725285;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 6://Daubiches6
wlen=12;
h[0]=0.111540743350;
h[1]=0.494623890398;
h[2]=0.751133908021;
h[3]=0.315250351709;
h[4]=-0.226264693965;
h[5]=-0.129766867567;
h[6]=0.097501605587;
h[7]=0.027522865530;
h[8]=-0.031582039318;
h[9]=0.000553842201;
h[10]=0.004777257511;
h[11]=-0.001077301085;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 7://Daubiches7
wlen=14;
h[0]=0.077852054085;
h[1]=0.396539319482;
h[2]=0.729132090846;
h[3]=0.469782287405;
h[4]=-0.143906003929;
h[5]=-0.224036184994;
h[6]=0.071309219267;
h[7]=0.080612609151;
h[8]=-0.038029936935;
h[9]=-0.016574541631;
h[10]=0.012550998556;
h[11]=0.000429577973;
h[12]=-0.001801640704;
h[13]=0.000353713800;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 8://Daubiches8
wlen=16;
h[0]=0.054415842243;
h[1]=0.312871590914;
h[2]=0.675630736297;
h[3]=0.585354683654;
h[4]=-0.015829105265;
h[5]=-0.284015542962;
h[6]=0.000472484574;
h[7]=0.128747426620;
h[8]=-0.017369301002;
h[9]=-0.044088253931;
h[10]=0.013981027917;
h[11]=0.008746094047;
h[12]=-0.004870352993;
h[13]=0.000391740373;
h[14]=0.000675449406;
h[15]=-0.000117476784;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 9://Daubiches9
wlen=18;
h[0]=0.038077947364;
h[1]=0.243834674613;
h[2]=0.604823123690;
h[3]=0.657288078051;
h[4]=0.133197385825;
h[5]=-0.293273783279;
h[6]=-0.096840783223;
h[7]=0.148540749338;
h[8]=0.030725681479;
h[9]=-0.067632829061;
h[10]=0.000250947115;
h[11]=0.022361662124;
h[12]=-0.004723204758;
h[13]=-0.004281503682;
h[14]=0.001847646883;
h[15]=0.000230385764;
h[16]=-0.000251963189;
h[17]=0.000039347320;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
case 10://db10
wlen=20;
h[0]=0.026670057901;
h[1]=0.188176800078;
h[2]=0.527201188932;
h[3]=0.688459039454;
h[4]=0.281172343661;
h[5]=-0.249846424327;
h[6]=-0.195946274377;
h[7]=0.127369340336;
h[8]=0.093057364604;
h[9]=-0.071394147166;
h[10]=-0.029457536822;
h[11]=0.033212674059;
h[12]=0.003606553567;
h[13]=-0.010733175483;
h[14]=0.001395351747;
h[15]=0.001992405295;
h[16]=-0.000685856695;
h[17]=-0.000116466855;
h[18]=0.000093588670;
h[19]=-0.000013264203;
for(i=0;i<wlen;i++)
g[i]=pow(-1,i)*h[-i+wlen-1];
break;
}
j=N;
flag[0]=0;
for(i=0;i<nscales;i++)//i的上限？？？？？
{
flag[i+1]=flag[i]+j;
sca[i]=j;
j=j/2;
}
switch(inv)
{
case 0:
mallat(a,d,h,g,wlen,nscales);
break;
case 1:
imallat(a,d,h,g,wlen,nscales);
break;
}
}
////////////////维格纳分布//////////////////////////////////////////
// 离散伪维格纳分布
// x-输入序列，输出为维格纳分布值，长度为(n2-n1+1)*nfft；
// N-输入序列x的长度；
// n1-维格纳分布的起始时刻，n1>=(L-1)/2；
// n2-维格纳分布的终止时刻，n2<N-(L-1)/2；
// L-分析窗的奇数宽度；
// nfft-FFT的长度，nfft>L;
// windowflag-窗函数类型;
extern "C" __declspec(dllexport)void WVD(float *x,int N,int n1,int n2,int L,int nfft,int windowflag)
{ float *xi,*real,*image, *gr, *gi,*y,*yi;
int i,j,k,B;
HANDLE hFloatxi,hFloatreal,hFloatimage,hFloatgr,hFloatgi,hFloaty,hFloatyi;
hFloatxi=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
xi=(float *)GlobalLock(hFloatxi);
hFloatreal=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,L*sizeof(float));
real=(float *)GlobalLock(hFloatreal);
hFloatimage=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,L*sizeof(float));
image=(float *)GlobalLock(hFloatimage);
hFloatgr=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,nfft*sizeof(float));
gr=(float *)GlobalLock(hFloatgr);
hFloatgi=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,nfft*sizeof(float));
gi=(float *)GlobalLock(hFloatgi);
hFloaty=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
y=(float *)GlobalLock(hFloaty);
hFloatyi=(HANDLE)GlobalAlloc(GMEM_MOVEABLE,N*sizeof(float));
yi=(float *)GlobalLock(hFloatyi);
for(i=0;i<N;i++)
{
xi[i]=0.0f;
}
Hilbert(x,xi,N);
for(i=0;i<N;i++)
{
y[i]=x[i];
yi[i]=xi[i];
}
for(i=n1;i<=n2;i++)
{
for(j=0;j<L;j++)
{ real[j]=y[i-(L-1)/2+j];
image[j]=yi[i-(L-1)/2+j];
}
windowed(real,L,windowflag);
windowed(image,L,windowflag);
for(k=0;k<=(L-1)/2;k++)
{
gr[k]=real[(L-1)/2+k]*real[(L-1)/2-k]+image[(L-1)/2+k]*image[(L-1)/2-k];
gi[k]=-real[(L-1)/2+k]*image[(L-1)/2-k]+image[(L-1)/2+k]*real[(L-1)/2-k];
}
for(k=(L-1)/2+1;k<nfft-(L-1)/2;k++)
gr[k]=gi[k]=0.0f;
for(k=1;k<=(L-1)/2;k++)
{
gr[nfft-k]=gr[k];
gi[nfft-k]=gi[k];
}
B=CalculatePower(nfft);
fft(gr,gi,B,0);
for(j=0;j<(nfft/2+1);j++)
{
x[(i-n1)*(nfft/2+1)+j]=gr[j];
// x[(i-n1)*(nfft/2+1)+j]=2*sqrt(gr[j]*gr[j]+gi[j]*gi[j]);;
}
}
GlobalUnlock(hFloatxi);
GlobalUnlock(hFloatreal);
GlobalUnlock(hFloatimage);
GlobalUnlock(hFloatgr);
GlobalUnlock(hFloatgi);
GlobalUnlock(hFloaty);
GlobalUnlock(hFloatyi);
GlobalFree(hFloatxi);
GlobalFree(hFloatreal);
GlobalFree(hFloatimage);
GlobalFree(hFloatgr);
GlobalFree(hFloatgi);
GlobalFree(hFloaty);
GlobalFree(hFloatyi);
}
///////////////////////巴氏滤波器/////////////////////////////////////
float FilterA[NODS/2],FilterB[NODS/2],FilterC[NODS/2];
float FilterD[NODS/2],FilterE[NODS/2];
/***************************************************/
/* Batterworth Digital Filter */
/***************************************************/
// fc:cut frequence fs:samplefrequence Nods:order of the filter
void LowPassFilter(float fc,float fs,int Nods)//b0=FilterA,b1=2FilterA,b2=FilterA;
{ float t,wcp,x,cs; //a0=1,a1=FilterB,a2=FilterC; smj
int k;
t=(float)(1.0/fs);
//wcp=(float)sin(fc*PI*t)/cos(fc*PI*t);
wcp=(float)tan(fc*PI*t);
for(k=1;k<=Nods/2;k++)
{
cs=(float)cos((2.*k+Nods-1.)*PI/2./Nods);
x=(float)(1./(1.+wcp*wcp-2.*wcp*cs));
FilterA[k]=wcp*wcp*x;
FilterB[k]=(float)(2.*(wcp*wcp-1.)*x);
FilterC[k]=(float)((1.+wcp*wcp+2.*wcp*cs)*x);
}
}
/***************************************************/
void HighPassFilter(float fc,float fs,int Nods)
{ float t,wcp,cs;
int k;
t=(float)(1.0/fs);
wcp=(float)tan(fc*PI*t);
for(k=1;k<=Nods/2;k++)
{
cs=(float)cos((2.*k+Nods-1.)*PI/2./Nods);
FilterA[k]=(float)(1./(1.+wcp*wcp-2.*wcp*cs));
FilterB[k]=(float)(2.*(wcp*wcp-1.)*FilterA[k]);
FilterC[k]=(float)((1.+wcp*wcp+2.*wcp*cs)*FilterA[k]);
}
}
/*********************************************************/
void BandPassFilter(float f1,float f2,float fs,int Nods)
{ float t,x,p,q,r,s,w1,w2,wc,cs;
int k;
t=(float)(1./fs);
w1=(float)tan(f1*PI*t);
w2=(float)tan(f2*PI*t);
wc=w2-w1;
q=(float)(wc*wc+2.*w1*w2);
s=w1*w1*w2*w2;
for(k=1;k<=Nods/2;k++)
{
cs=(float)cos((2.*k+Nods-1.)*PI/2./Nods);
p=(float)(-2.*wc*cs);
r=p*w1*w2;
x=1+p+q+r+s;
FilterA[k]=wc*wc/x;
FilterB[k]=(-4-2*p+2*r+4*s)/x;
FilterC[k]=(6-2*q+6*s)/x;
FilterD[k]=(-4+2*p-2*r+4*s)/x;
FilterE[k]=(float)((1.-p+q-r+s)/x);
}
}
/*********************************************************/
// 低/高通滤波
// x--输入序列，输出为滤波后结果；
// N--序列长度；
// Fc--截止频率；
// Fs--采样频率；
// Nods--ButterWorth滤波器的阶数；
// LorH--滤波选择，-1为低通，其它为高通；
extern "C" __declspec(dllexport)void LorH_Pass(float *x,int N,float fc,float fs,int Nods,int LorH)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int k,i;
float *y;
float cc;
y=(float *)malloc(N*sizeof(float));
if(LorH==-1)
{ LowPassFilter(fc,fs,Nods);
cc=2.;
}
else
{ HighPassFilter(fc,fs,Nods);
cc=-2.;
}
for(k=1;k<=Nods/2;k++)
{
y[0]=FilterA[k]*x[0];
y[1]=FilterA[k]*(x[1]+cc*x[0])-FilterB[k]*y[0];
y[2]=FilterA[k]*(x[2]+cc*x[1])-FilterB[k]*y[1]-FilterC[k]*y[0];
for(i=2;i<N;i++)
y[i]=FilterA[k]*(x[i]+cc*x[i-1]+x[i-2])-FilterB[k]*y[i-1]-FilterC[k]*y[i-2];
for(i=0;i<N;i++)
x[i]=y[i];
}
free(y);
}
/*********************************************************/
// 带通滤波
// x--输入序列，输出为滤波后结果；
// N--序列长度；
// fLeft--下截止频率；
// fRight--上截止频率；
// fs--采样频率；
// Nods--ButterWorth滤波器的阶数；
extern "C" __declspec(dllexport)void BandPass(float *x,int N,float fLeft,float fRight,float fs,int Nods)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int k,i;
float *y;
BandPassFilter(fLeft,fRight,fs,Nods);
y=(float *)malloc(N*sizeof(float));
for(k=1;k<=Nods/2;k++)
{
y[0]=FilterA[k]*x[0];
y[1]=FilterA[k]*x[1]-FilterB[k]*y[0];
y[2]=FilterA[k]*(x[2]-2*x[0])-FilterB[k]*y[1]-FilterC[k]*y[0];
y[3]=FilterA[k]*(x[3]-2*x[1])-FilterB[k]*y[2]-FilterC[k]*y[1]-FilterD[k]*y[0];
for(i=4;i<N;i++)
y[i]=FilterA[k]*(x[i]-2*x[i-2]+x[i-4])-FilterB[k]*y[i-1]-FilterC[k]*y[i-2]-FilterD[k]*y[i-3]-FilterE[k]*y[i-4];
for(i=0;i<N;i++)
x[i]=y[i];
}
free(y);
}
// 带阻滤波
// x--输入序列，输出为滤波后结果；
// N--序列长度；
// fLeft--下截止频率；
// fRight--上截止频率；
// fs--采样频率；
// Nods--ButterWorth滤波器的阶数；
extern "C" __declspec(dllexport)void BandStop(float *x,int N,float fLeft,float fRight,float fs,int Nods)
{
AFX_MANAGE_STATE(AfxGetStaticModuleState());
int k,i;
float t,w1,w2,wc;
float b0,b1,b2,b3,b4;
float *y;
BandPassFilter(fLeft,fRight,fs,Nods);
t=(float)(1./fs);
w1=(float)tan(fLeft*PI*t);
w2=(float)tan(fRight*PI*t);
wc=(w2-w1)*(w2-w1);
b0=(float)(1.+2.*w1*w2+w1*w1*w2*w2);
b1=(float)(-4.+4.*w1*w1*w2*w2);
b2=(float)(6.-4.*w1*w2+6.*w1*w1*w2*w2);
b3=b1;
b4=b0;
y=(float *)malloc(N*sizeof(float));
for(k=1;k<=Nods/2;k++)
{
y[0]=FilterA[k]*b0*x[0]/wc;
y[1]=FilterA[k]*(b0*x[1]+b1*x[0])/wc-FilterB[k]*y[0];
y[2]=FilterA[k]*(b0*x[2]+b1*x[1]+b2*x[0])/wc-FilterB[k]*y[1]-FilterC[k]*y[0];
y[3]=FilterA[k]*(b0*x[3]+b1*x[2]+b2*x[1]+b3*x[0])/wc-FilterB[k]*y[2]-FilterC[k]*y[1]-FilterD[k]*y[0];
for(i=4;i<N;i++)
y[i]=FilterA[k]*(b0*x[i]+b1*x[i-1]+b2*x[i-2]+b3*x[i-3]+b4*x[i-4])/wc-FilterB[k]*y[i-1]-FilterC[k]*y[i-2]-FilterD[k]*y[i-3]-FilterE[k]*y[i-4];
for(i=0;i<N;i++)
x[i]=y[i];
}
free(y);
}
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?