边缘滤波问题c语言,关于 filtfilt 函数（c语言版）边缘效应的消除问题

已结贴√

问题点数：20 回复次数：1

关于 filtfilt 函数(c语言版)边缘效应的消除问题

//filtfilt函数

int filtfilt(double* x, double* y, int xlen, double* a, double* b, int nfilt)

{

int nfact;

int tlen;    //length of tx

int i;

double *tx,*tx1,*p,*t,*end;

double *sp,*tvec,*zi;

double tmp,tmp1;

nfact=nfilt-1;    //3*nfact: length of edge transients

if(xlen<=3*nfact || nfilt<2) return -1;

tlen=6*nfact+xlen;

tx=(double *)malloc(tlen*sizeof(double));

tx1=(double *)malloc(tlen*sizeof(double));

sp=(double *)malloc( sizeof(double) * nfact * nfact );

tvec=(double *)malloc( sizeof(double) * nfact );

zi=(double *)malloc( sizeof(double) * nfact );

if( !tx || !tx1 || !sp || !tvec || !zi )

{

free(tx);

free(tx1);

free(sp);

free(tvec);

free(zi);

return 1;

}

tmp=x[0];

for(p=x+3*nfact,t=tx;p>x;--p,++t) *t=2*tmp-*p;

for(end=x+xlen;p

tmp=x[xlen-1];

for(end=tx+tlen,p-=2;t

//now tx is ok.

end = sp + nfact*nfact;

p=sp;

while(p

sp[0]=1.0+a[1];

for(i=1;i

sp[i*nfact]=a[i+1];

sp[i*nfact+i]=1.0L;

sp[(i-1)*nfact+i]=-1.0L;

}

for(i=0;i

tvec[i]=b[i+1]-a[i+1]*b[0];

}

if(rinv(sp,nfact)){

free(zi);

zi=NULL;

}

else{

trmul(sp,tvec,zi,nfact,nfact,1);

}//zi is ok

free(sp);free(tvec);

//filtering tx, save it in tx1

tmp1=tx[0];

if(zi)

for( p=zi,end=zi+nfact; p

filter(tx,tx1,tlen,a,b,nfilt,zi);

//reverse tx1

for( p=tx1,end=tx1+tlen-1; p

tmp = *p;

*p = *end;

*end = tmp;

}

//filter again

tmp1 = (*tx1)/tmp1;

if(zi)

for( p=zi,end=zi+nfact; p

filter(tx1,tx,tlen,a,b,nfilt,zi);

//reverse to y

end = y+xlen;

p = tx+3*nfact+xlen-1;

while(y

*y++ = *p--;

}

free(zi);

free(tx);

free(tx1);

return 0;

}

这段代码是c语言实现零相移滤波的滤波函数，但是验证的结果最终波形开头和结尾有很大的波动，而中间的跟matlab验证是重合的，里面有对开头和结尾数据进行延拓，但是我不明白原理，有人能解答一下吗