贡献一个由滤波器系数计算滤波后数据的程序
看到论坛上有人问由滤波器系数求滤波后数据的问题（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;
double *xdelayed=new double;
for(int i=0;i<len_x;i++)
{
v=coeff_b*x;
xdelayed=0;filter_x=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=x;
}
for(int j=0;j<len_x;j++)
{
v=v+coeff_b*xdelayed;
}
xdelayed=0;
}
}
if(len_a>1)
{
double *ac=new double;
for(int i=0;i<len_a-1;i++)
ac=-1*coeff_a;
for(int i1=0;i1<len_x;i1++)
{
double t=v;
for(int j=0;j<len_a-1;j++)
{
if(i1>j)
{
t=t+ac*filter_x;
}
}
filter_x=t;
}
}
else
for(int i=0;i<len_x;i++)
{
filter_x=v;
}
}
顺便把matlab源程序也贴一下：
function = 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 = ;
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
一点小改动，防止内存泄露。
在函数结尾加上
delete []v;
delete []xdelayed;
在if(len_a>1)
{
}结尾处也加上
delete []ac;
页:
[1]