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];
y1 = filterimp(a, b, 0, 4, x, 1);

figure('name', '平均值滤波器'); grid on; hold on;
plot(t,x ,'r')
plot(t,y ,'b*')
plot(t,y1 ,'g')
legend('Input Data','Filter','Filter Imp');

测试示例1，矩阵滤波测试

%对矩阵行进行滤波
figure('name', '对矩阵行进行滤波'); grid on; hold on;
x = rand(2,15);%2*15的矩陣
b = 1;
a = [1 -0.2];

y = filter(b,a,x,[],2);

t = 0:length(x)-1;  %index vector

plot(t,x(1,:),'r')
hold on
plot(t,y(1,:),'b*')
plot(t,y(2,:),'m--o')

y1 = filterimp(a, b, 1, 0, x,2);
plot(t,y1(1,:) ,'g')
plot(t,y1(2,:) ,'k')
legend('Input Data','Filter1','Filter2' ,'Filterd Imp1','Filter Imp2');

测试示例3: a b 为滤波数组

%a b 滤波
b = [2,3];
a = [1,0.2];
x = [1,2,3,4,5,6,7,8,1,2,3,4];

y = filter(b,a,x);
y1 = filterimp(a, b, 1, 1, x,1);
figure('name', 'a b 滤波'); grid on; hold on;
plot(x ,'r')
plot(y ,'b*')
plot(y1 ,'g')
legend('Input Data','Filter','Filter Imp');