1、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 describedby vectors A and B to create t
2、he 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
3、 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 t
4、imesthe filter order, defined as max(length(b)-1,length(a)-1).Note that FILTFILT should not be used with differentiator and Hilbert FIRfilters, since the operation of these filters depends heavily on theirphase response.See also filter.Reference page in Help browserdoc filtfilt零相移滤波 filtfilt，操作上不复杂，
5、原理上小小有点复杂。它主要是先找到一个“合理的”初始状态 Zi，使得无论先从反向滤波还是先从正向滤波结果一致，然后 filter 一次，逆向，再filter 一次，再逆向，就是结果了。这里面包括 3 个要点：1. Zi 的确定。2. 边缘效应的消除。3. 正反向滤波的数学原理。对于要点 1，可以参阅 Fredrik Gustafsson 的论文 Determining the Initial States in Forward-backward Filtering 的数学证明。要点 2，filtfilt 对数据两头做了镜像拓延，最后滤完波再截掉头尾。要点 3，这个貌似不大难推#include
6、#include “filter.h“#include “mulMat.h“#include “invMat.h“intfiltfilt(const double* x, double* y, int xlen, double* a, double* b, int nfilt)int nfact;int tlen; /length of txint i;double *tx,*tx1,*p,*t,*end;double *sp,*tvec,*zi;double tmp,tmp1;nfact=nfilt-1; /3*nfact: length of edge transientsif(xlenx
7、;-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 spsp0=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
8、-ai+1*b0;if(invMat(sp,nfact)free(zi);zi=NULL;elsemulMat(sp,tvec,zi,nfact,nfact,1);/zi is okfree(sp);free(tvec);/filtering tx, save it in tx1tmp1=tx0;if(zi)for( p=zi,end=zi+nfact; pend;) *(p+) *= tmp1;filter(tx,tx1,tlen,a,b,nfilt,zi);/reverse tx1for( p=tx1,end=tx1+tlen-1; pend; p+,end-)tmp = *p;*p = *end;*end = tmp;/filter againtmp1 = (*tx1)/tmp1;if(zi)for( p=zi,end=zi+nfact; pend;) *(p+) *= tmp1;filter(tx1,tx,tlen,a,b,nfilt,zi);/reverse to yend = y+xlen;p = tx+3*nfact+xlen-1;while(yend)*y+ = *p-;free(zi);free(tx);free(tx1);return 0;
