matlab filter函数c实现

 x输入，y输出，xlen长度，a为分子系数，b分母系数，nfilt滤波器阶数

#include <stdio.h>
#include <math.h>
//#include <malloc.h>
#include <string.h>
//#include "filtfilt.h"
#define EPS 0.000001
//filter函数
void filter(const double* x, double* y, int xlen, double* a, double* b, int nfilt)  
{       
    double tmp;  
    int i,j;  
  
    //normalization  
    if( (*a-1.0>EPS) || (*a-1.0<-EPS) )  
    {  
        tmp=*a;  
        for(i=0;i<nfilt;i++)  
        {  
            b[i]/=tmp;  
            a[i]/=tmp;  
        }  
    }  
  
    memset(y,0,xlen*sizeof(double));  
  
    a[0]=0.0;  
    for(i=0;i<xlen;i++)  
    {  
        for(j=0;i>=j&&j<nfilt;j++)  
        {  
            y[i] += (b[j]*x[i-j]-a[j]*y[i-j]);  
        }  
       // if(zi&&i<nfilt-1) y[i] += zi[i];   //zi[i],查分方程可以不计，滤波的表达式 1.[y, zf] = filter(b ,a, X) 2.[y, zf] = filter(b, a, X, zi)  
    }  
    a[0]=1.0;  
  
}