贡献一个由滤波器系数计算滤波后数据的程序

impulse

看到论坛上有人问由滤波器系数求滤波后数据的问题（matlab中有相应的内建函数filter），目前只能针对一组滤波器系数b，不能针对两组滤波器系数b和a。（我也是从网上搜到的matlab程序，把它改变为C++）。
void filter(double *x, int len_x, double *coeff_b, int len_b, double *coeff_a, int len_a,double *filter_x)
{
        double *v=new double[len_x];
        double *xdelayed=new double[len_x];
        for(int i=0;i<len_x;i++)
        {
                v[i]=coeff_b[0]*x[i];
                xdelayed[i]=0;filter_x[i]=0;
        }
        if(len_b>1)
        {
                int mi=min(len_x,len_b);
                for(int i=1;i<mi;i++)
                {
                        for(int k=0;k<len_x-i;k++)
                        {
                                xdelayed[i+k]=x[k];                                       
                        }
                        for(int j=0;j<len_x;j++)
                        {
                                v[j]=v[j]+coeff_b[i]*xdelayed[j];
                               
                        }
                        xdelayed[i]=0;
                }
        }
        if(len_a>1)
        {
                double *ac=new double[len_a-1];
                for(int i=0;i<len_a-1;i++)
                        ac[i]=-1*coeff_a[i+1];
               
                for(int i1=0;i1<len_x;i1++)
                {
                        double t=v[i1];
                        for(int j=0;j<len_a-1;j++)
                        {
                                if(i1>j)
                                {
                                        t=t+ac[j]*filter_x[i1-j-1];
                                }
                        }
                        filter_x[i1]=t;
                }
        }
        else
                for(int i=0;i<len_x;i++)
                {
                        filter_x[i]=v[i];
                }       
               
}
顺便把matlab源程序也贴一下：
function [y] = filterslow(B,A,x)
% FILTERSLOW: Filter x to produce y = (B/A) x .
%       Equivalent to 'y = filter(B,A,x)' using
%       a slow (but tutorial) method.

NB = length(B);
NA = length(A);
Nx = length(x);

xv = x(:); % ensure column vector

% do the FIR part using vector processing:
v = B(1)*xv;
if NB>1
  for i=2:min(NB,Nx)
    xdelayed = [zeros(i-1,1); xv(1:Nx-i+1)];
    v = v + B(i)*xdelayed;
  end;
end; % fir part done, sitting in v

% The feedback part is intrinsically scalar,
% so this loop is where we spend a lot of time.
y = zeros(length(x),1); % pre-allocate y
ac = - A(2:NA);
for i=1:Nx, % loop over input samples
  t=v(i);   % initialize accumulator
  if NA>1,
    for j=1:NA-1
      if i>j,
        t=t+ac(j)*y(i-j)
      %else y(i-j) = 0
      end;
    end;
  end;
  y(i)=t;
  clc;
end;

y = reshape(y,size(x)); % in case x was a row vector

impulse

一点小改动，防止内存泄露。
在函数结尾加上
delete []v;
delete []xdelayed;
在if(len_a>1)
{

}结尾处也加上
delete []ac;