注册 登录
编程论坛 数据结构与算法

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

fu563048951 发布于 2014-07-18 14:44, 1327 次点击
//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验证是重合的,里面有对开头和结尾数据进行延拓,但是我不明白原理,有人能解答一下吗
1 回复
#2
yuccn2014-07-24 12:35
同样不明白
1