filter 滤波代码实现

filter滤波函数实现

%用代码实现filter
%na = length(a) - 1
%nb = length(b) -1
%xn 为矩阵x的维数
function res = filterimp(a ,b ,na ,nb ,x ,xn)
% a(1)*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)
y = [];

for row = 1:xn
    y(row,1) = (b(1) * x(row,1))/a(1);
    for n = 2 : length(x)
        z1 = 0;
        for j = 1:(nb + 1)
            if j <= n
                z1 = z1 + b(j)*x(row,n-j+1);
            end
        end
        z2 = 0;
        for i = 2:(na+1)
            if i <= n
                z2 = a(i)*y(row,n-i+1);
            end
        end
        y(row,n) = z1 - z2;
        y(row,n) = y(row,n)/a(1);
    end
end

res = y;

测试示例1，滤波正弦波噪声

t = linspace(-pi,pi,100);
rng default  %initialize random number generator
x = sin(t) + 0.25*rand(size(t));

windowSize = 5; 
b = (1/windowSize)*ones(1,windowSize);%[0.2 0.2 0.2 0.2 0.2]
a = 1;
y = filter(b,a,x);

a = [1