EEPW论坛

最近在做轴承故障检测，C语言编程，在滤波部分遇到了问题，由于本人对滤波知之甚少，想求助：

这是测得的振动信号，在这个基础上我分别用C语言与LabVIEW对其做了滤波， 都是用iir巴特沃兹滤波， 频域上两者结果一样（我就不放图了），但时域差距很大。

这是LabVIEW做得滤波之后的时域图像看上去很正常，iir巴特沃兹滤波3阶

这是C语言同样采用iir巴特沃兹滤波的时域结果，非常不像。为什么会出现这样的问题。C语言是用的双线性变换

void iirbcf(short ifilt,short band,short ns,short n,double f1,double f2,double f3,double f4,double db,double b[],double a[]) { short k; double omega,lamda,epslon,fl,fh; double d[5],c[5]; if((band==1)||(band==4)) fl=f1; if((band==2)||(band==3)) fl=f2; if(band<=3) fh=f3; if(band==4) fh=f4; if(ifilt<3) { switch (band) { case 1: case 2: { omega=warp(f2)/warp(f1); break; } case 3: { omega=omin(bpsub(warp(f1),fh,fl),bpsub(warp(f4),fh,fl)); break; } case 4: { omega=omin(1.0/bpsub(warp(f2),fh,fl),1.0/bpsub(warp(f3),fh,fl));} } lamda=pow(10.0,(db/20.0)); epslon=lamda/cosh(2*ns*cosh1(omega)); } for(k=0;k=0;i--) { if((c[i]!=0.0)||(d[i]!=0.0)) break; } m=i; switch (band) { case 1: case 2: { n2=m; n1=n2+1; if(band==2) { for(i=0;i<=m/2;i++) { tmp=d[i]; d[i]=d[m-i]; d[m-i]=tmp; tmp=c[i]; c[i]=c[m-i]; c[m-i]=tmp; } } for(i=0;i<=m;i++) { d[i]=d[i]/pow(w1,i); c[i]=c[i]/pow(w1,i); } break; } case 3: case 4: { n2=2*m; n1=n2+1; work=malloc(n1*n1*sizeof(double)); w2=tan(pi*fhn); w=w2-w1; w0=w1*w2; if(band==4) { for(i=0;i<=m/2;i++) { tmp=d[i]; d[i]=d[m-i]; d[m-i]=tmp; tmp=c[i]; c[i]=c[m-i]; c[m-i]=tmp; } } for(i=0;i<=n2;i++) { work[0*n1+i]=0.0; work[1*n1+i]=0.0; } for(i=0;i<=m;i++) { tmpd=d[i]*pow(w,(m-i)); tmpc=c[i]*pow(w,(m-i)); for(k=0;k<=i;k++) { ls=m+i-2*k; tmp=combin(i,i)/(combin(k,k)*combin(i-k,i-k)); work[0*n1+ls]+=tmpd*pow(w0,k)*tmp; work[1*n1+ls]+=tmpc*pow(w0,k)*tmp; } } for(i=0;i<=n2;i++) { d[i]=work[0*n1+i]; c[i]=work[1*n1+i]; } free(work); } } bilinear(d,c,b,a,n); } double combin(short i1,short i2) { short i; double s; s=1.0; if(i2==0) return (s); for(i=i1;i>(i1-i2);i--) { s*=i; } return (s); } void bilinear(double d[],double c[],double b[],double a[],short n) { short i,j,n1; double sum,atmp,scale, *temp; n1=n+1; temp=malloc(n1*n1*sizeof(double)); for(j=0;j<=n;j++) { temp[j*n1+0]=1.0; } sum=1.0; for(i=1;i<=n;i++) { sum=sum*(double)(n-i+1)/(double)i; temp[0*n1+i]=sum; } for(i=1;i<=n;i++) for(j=1;j<=n;j++) { temp[j*n1+i]=temp[(j-1)*n1+i]-temp[j*n1+i-1]-temp[(j-1)*n1+i-1]; } for(i=n;i>=0;i--) { b[i]=0.0; atmp=0.0; for(j=0;j<=n;j++) { b[i]=b[i]+temp[j*n1+i]*d[j]; atmp=atmp+temp[j*n1+i]*c[j]; } scale=atmp; if(i!=0) a[i]=atmp; } for(i=0;i<=n;i++) { b[i]=b[i]/scale; a[i]=a[i]/scale; } a[0]=1.0; free(temp); }

这是C语言滤波程序，我只用了里面的巴特沃兹