Sample-Level Implementation in Matlab

For completeness, a direct matlab implementation of the built-in filter function (Eq.(3.3)) is given in Fig.3.2. While this code is useful for study, it is far slower than the built-in filter function. As a specific example, filtering samples of data using an order 100 filter on a 900MHz Athlon PC required 0.01 seconds for filter and 10.4 seconds for filterslow. Thus, filter was over a thousand times faster than filterslow in this case. The complete test is given in the following matlab listing:

x = rand(10000,1); % random input signal
B = rand(101,1);   % random coefficients
A = [1;0.001*rand(100,1)]; % random but probably stable
tic; yf=filter(B,A,x); ft=toc
tic; yfs=filterslow(B,A,x); fst=toc

1. recursive feedback cannot be ``vectorized'' in general, and
2. built-in functions such as filter are written in C, precompiled, and linked with the main program.

function [y] = 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 = [zeros(i-1,1); xv(1:Nx-i+1)];
    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; 
end; 

y = reshape(y,size(x)); % in case x was a row vector

Order a Hardcopy of Introduction to Digital Filters