请教滤波问题

angle3839

小弟用matlab的4阶巴特沃斯滤波器滤波：
滤波系数：B=[4.383200469249005e-010    1.753280187699602e-009   2.629920281549403e-009    1.753280187699602e-009   4.383200469249005e-010 ];
A =[1.000000000000000e+000   -3.976015181233625e+000   5.928332760985220e+000   -3.928617967640213e+000  9.763003949017387e-001 ];
滤波函数是这么写的：
%myfilter
load 'y2lf_d.txt'
x=zeros(length(y2lf_d),1);
y=x;
x=y2lf_d;
B=[4.383200469249005e-010    1.753280187699602e-009   2.629920281549403e-009    1.753280187699602e-009   4.383200469249005e-010 ];
A =[1.000000000000000e+000   -3.976015181233625e+000   5.928332760985220e+000   -3.928617967640213e+000  9.763003949017387e-001 ];
N=length(B);
for m=1:length(y2lf_d);
switch m
    case 1:N-1
        y(1)=B(1)*x(1),format long e;
        y(2)=B(1)*x(2)+B(2)*x(1)-A(2)*y(1),format long e;
        y(3)=B(1)*x(3)+B(2)*x(2)+B(3)*x(1)-(A(2)*y(2)+A(3)*y(1)),format long e;
        y(4)=B(1)*x(4)+B(2)*x(3)+B(3)*x(2)+B(4)*x(1)-(A(2)*y(3)+A(3)*y(2)+A(4)*y(1)),format long e;
        %y(5)=B(1)*x(5)+B(2)*x(4)+B(3)*x(3)+B(4)*x(2)+B(5)*x(1)-(A(2)*y(4)+A(3)*y(3)+A(4)*y(2)+A(5)*y(1)),format long e;
    case N+1:length(y2lf_d)
        y(m)=B(1)*x(m)+B(2)*x(m-1)+B(3)*x(m-2)+B(4)*x(m-3)+B(5)*x(m-4)-(A(2)*y(m-1)+A(3)*y(m-2)+A(4)*y(m-3)+A(5)*y(m-4));
    otherwise
        
end
end
，滤波全成零了，如下图：

                               
登录/注册后可看大图

而用它自带的filter函数则滤波效果挺好，如下图：

                               
登录/注册后可看大图

                               
登录/注册后可看大图

为什么？
原理不是一样的吗？我那个就是按照滤波的定义写的呀。
哪位帮帮忙，谢谢啊！

[ 本帖最后由 xinyuxf 于 2007-4-30 16:35 编辑 ]

songzy41

存在的错误有：
1，对y(1)~y(4)的计算，不应在for内；
2，笫2个循环从N开始。
滤波过程可直接写成：
        y(1)=B(1)*x(1),format long e;
        y(2)=B(1)*x(2)+B(2)*x(1)-A(2)*y(1),format long e;
        y(3)=B(1)*x(3)+B(2)*x(2)+B(3)*x(1)-(A(2)*y(2)+A(3)*y(1)),format long e;
        y(4)=B(1)*x(4)+B(2)*x(3)+B(3)*x(2)+B(4)*x(1)-(A(2)*y(3)+A(3)*y(2)+A(4)*y(1)),format long e;      
for m=N:length(y2lf_d)
        y(m)=B(1)*x(m)+B(2)*x(m-1)+B(3)*x(m-2)+B(4)*x(m-3)+B(5)*x(m-4)-(A(2)*y(m-1)+A(3)*y(m-2)+A(4)*y(m-3)+A(5)*y(m-4));   
end

angle3839

谢谢，谢谢！明白了。

angle3839

怎么评分，我也想评，怎么办？是不是我没有权限？

eight

原帖由 angle3839 于 2007-4-30 18:38 发表

                               
登录/注册后可看大图

怎么评分，我也想评，怎么办？是不是我没有权限？

请查看导航条中的“我的权限”以及论坛关于普通会员等级的公告