filtfilt函数的c语言实现.doc

y,zf = filter(b,a,X)y,zf = filter(b,a,X,zi)其中b和a就是差分方程的系数，X是输入向量，zi是“初始状态”。可能这么说明还是不很清晰，那么请看图（注意，a(1)为1，这个可以通过差分方程所有系数除以a(1)所得）：FILTFILT Zero-phase forward and reverse digital filtering. Y = FILTFILT(B, A, X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y. The filter is described by the difference equation: y(n) = b(1)*x(n) + b(2)*x(n-1) + . + b(nb+1)*x(n-nb) - a(2)*y(n-1) - . - a(na+1)*y(n-na) After filtering in the forward direction, the filtered sequence is then reversed and run back through the filter; Y is the time reverse of the output of the second filtering operation. The result has precisely zero phase distortion and magnitude modified by the square of the filters magnitude response. Care is taken to minimize startup and ending transients by matching initial conditions. The length of the input x must be more than three times the filter order, defined as max(length(b)-1,length(a)-1). Note that FILTFILT should not be used with differentiator and Hilbert FIR filters, since the operation of these filters depends heavily on their phase response. See also filter. Reference page in Help browser doc filtfilt零相移滤波filtfilt，操作上不复杂，原理上小小有点复杂。它主要是先找到一个“合理的”初始状态Zi，使得无论先从反向滤波还是先从正向滤波结果一致，然后filter一次，逆向，再filter一次，再逆向，就是结果了。这里面包括3个要点：1. Zi的确定。2. 边缘效应的消除。3. 正反向滤波的数学原理。对于要点1，可以参阅Fredrik Gustafsson的论文Determining the Initial States in Forward-backward Filtering的数学证明。要点2，filtfilt对数据两头做了镜像拓延，最后滤完波再截掉头尾。要点3，这个貌似不大难推#include #include filter.h#include mulMat.h#include invMat.hintfiltfilt(const 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 | nfiltx;-p,+t) *t=2.0*tmp-*p; for(end=x+xlen;pend;+p,+t) *t=*p; tmp=xxlen-1; for(end=tx+tlen,p-=2;tend;-p,+t) *t=2.0*tmp-*p; /now tx is ok. end = sp + nfact*nfact; p=sp; while(pend) *p+ = 0.0L; /clear sp sp0=1.0+a1; for(i=1;infact;i+) spi*nfact=ai+1; spi*nfact+i=1.0L; sp(i-1)*nfact+i=-1.0L; for(i=0;infact;i+) tveci=bi+1-ai+1*b0; if(invMat(sp,nfact) free(zi); zi=NULL; else mulMat(sp,tvec,zi,nfact,nfact,1); /zi is ok free(sp);free(tvec); /filtering tx, save it in tx1 tmp1=tx0; if(zi) for( p=zi,end=zi+nfact; pend;) *(p+) *= tmp1; filter(tx,tx1,tlen,a,b,nfilt,zi); /reverse tx1 for( p=tx1,end=tx1+tlen-1; pend; p+,end-) tmp = *p; *p = *end; *end = tmp; /filter again tmp1 = (*tx1)/tmp1; if(zi) for( p=zi,end=zi+nfact; pend;) *(p+) *= tmp1; filter(tx1,tx,