| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 1327 人关注过本帖
标题:关于 filtfilt 函数(c语言版)边缘效应的消除问题
只看楼主 加入收藏
fu563048951
Rank: 2
等 级:论坛游民
帖 子:18
专家分:20
注 册:2012-2-19
结帖率:75%
收藏
已结贴  问题点数: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<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验证是重合的,里面有对开头和结尾数据进行延拓,但是我不明白原理,有人能解答一下吗
搜索更多相关主题的帖子: c语言 return double 
2014-07-18 14:44
yuccn
Rank: 16Rank: 16Rank: 16Rank: 16
来 自:何方
等 级:版主
威 望:167
帖 子:6815
专家分:42393
注 册:2010-12-16
收藏
得分:20 
同样不明白

我行我乐
公众号:逻辑客栈
我的博客:
https://blog.yuccn. net
2014-07-24 12:35
快速回复:关于 filtfilt 函数(c语言版)边缘效应的消除问题
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.025026 second(s), 9 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved