How can I write the Matlab "filter"-function myself?

Question

I would like to use a Butterworth filter on a 1D-Signal. In Matlab the script would look like this:

 f=100;
 f_cutoff = 20; 
 fnorm =f_cutoff/(f/2);
 [b,a] = butter(8,fnorm,'low');
 filteredData = filter(b,a,rawData); % I want to write this myself

Now I don't want to directly use the filter-function given in Matlab but write it myself. In the Matlab documentation it's described as follows:

The filter function is implemented as a direct form II transposed structure,

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)

where n-1 is the filter order, which handles both FIR and IIR filters [1], na is the feedback filter order, and nb is the feedforward filter order.

So I've already tried to write the function like that:

f=100;
f_cutoff = 20; 
fnorm =f_cutoff/(f/2);
[b,a] = butter(8,fnorm,'low');
for n = 9:size(rawData,1)
    filteredData(n,1) = b(1)*n + b(2)*(n-1) + b(3)*(n-2) + b(4)*(n-3) + b(5)*(n-4) ...
                      - a(2)*rawData(n-1,1) - a(3)*rawData(n-2,1) - a(4)*rawData(n-3,1) - a(5)*accel(n-4,1);
end

But that's not working. Can you please help me? What am I doing wrong?

Sincerely, Cerdo

PS: the filter documentation can be foud here: http://www.mathworks.de/de/help/matlab/ref/filter.html#f83-1015962 when expanding More About -> Algorithms

Answer 1

I have found a text described the Direct Form II Transposed used in the Matlab filter function and it works perfectly. See script below. Other implementations are also available but with error of around 1e-15, you'll see this by running the script yourself.

%% Specification of the Linear Chebysev filters
clc;clear all;close all
ord = 5; %System order (from 1 to 5)
[bq,aq] = cheby1(ord,2,0.2);theta = [bq aq(2:end)]';
figure;zplane(bq,aq); % Z-Pole/Zeros
u = [ones(40,1); zeros(40,1)];
%% Naive implementation of the basic algorithm
y0 = filter(bq,aq,u); % Built-in filter
b = fliplr(bq);a = fliplr(aq);a(end) = [];
y1 = zeros(40,1);pad = zeros (ord,1);
yp = [pad; y1(:)];up = [pad; u(:)];
for i = 1:length(u)
    yp(i+ord) = sum(b(:).*up(i:i+ord))-sum(a(:).*yp(i:i+ord-1));
end
y1 = yp(ord+1:end); % Naive implementation
err = y0(:)-y1(:);
figure
plot(y0,'r')
hold on
plot(y1,'*g')
xlabel('Time')
ylabel('Response')
legend('My code','Built-in filter')
figure
plot(err)
xlabel('Time')
ylabel('Error')
%% Direct Form II Transposed  
% Direct realization of rational transfer functions
% trps: 0 for direct realization, 1 for transposed realisation 
% b,a: Numerator and denominator
% x:   Input sequence
% y:   Output sequence
% u:   Internal states buffer

trps = 1;
b=theta(1:ord+1);
a=theta(ord+2:end);
y2=zeros(size(u));
x=zeros(ord,1);
%%
if trps==1
    for i=1:length(u)
        y2(i)=b(1)*u(i)+x(1);
        x=[x(2:ord);0];
        x=x+b(2:end)*u(i)-a*y2(i);
    end
else
    for i=1:length(u)
        xnew=u(i)-sum(x(1:ord).*a);
        x=[xnew,x];
        y2(i)=sum(x(1:ord+1).*b);
        x=x(1:ord);
    end
end
%%
err = y2 - filter(bq,aq,u);
figure
plot(y0,'r')
hold on
plot(y2,'*g')
xlabel('Time')
ylabel('Response')
legend('Form II Transposed','Built-in filter')
figure
plot(err)
xlabel('Time')
ylabel('Error')
% end