关于 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<end;++p,++t) *t=*p;
tmp=x[xlen-1];
for(end=tx+tlen,p-=2;t<end;--p,++t) *t=2*tmp-*p;
//now tx is ok.
end = sp + nfact*nfact;
p=sp;
while(p<end) *p++ = 0.0L; //clear sp
sp[0]=1.0+a[1];
for(i=1;i<nfact;i++){
sp[i*nfact]=a[i+1];
sp[i*nfact+i]=1.0L;
sp[(i-1)*nfact+i]=-1.0L;
}
for(i=0;i<nfact;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<end;) *(p++) *= tmp1;
filter(tx,tx1,tlen,a,b,nfilt,zi);
//reverse tx1
for( p=tx1,end=tx1+tlen-1; p<end; p++,end--){
tmp = *p;
*p = *end;
*end = tmp;
}
//filter again
tmp1 = (*tx1)/tmp1;
if(zi)
for( p=zi,end=zi+nfact; p<end;) *(p++) *= tmp1;
filter(tx1,tx,tlen,a,b,nfilt,zi);
//reverse to y
end = y+xlen;
p = tx+3*nfact+xlen-1;
while(y<end){
*y++ = *p--;
}
free(zi);
free(tx);
free(tx1);
return 0;
}
这段代码是c语言实现零相移滤波的滤波函数,但是验证的结果最终波形开头和结尾有很大的波动,而中间的跟matlab验证是重合的,里面有对开头和结尾数据进行延拓,但是我不明白原理,有人能解答一下吗