%-------------------------------------------------------------------------
%     *** 2.161 Signal Processing - Continuous and Discrete  ***
%                      Tutorial MATLAB Function
%
%  iirdf1 - Demonstration IIR Direct Form I filter implementation.
% 
% Source:   Class handout:  Direct Form Digital Filter Structures
%
% Usage :   1) Initialization:
%              y = iirdf1('initial', b, a)
%                where b, a are the numerator and denominator polynomial 
%                coefficients.  Example:
%                   [b,a] = butter(7,0.4);
%                    y = iirdf1('initial',b,a);
%                Note: iirdf1 returns y = 0 for initialization
%           2) Filtering:
%              y_out = iirdf1(f_{in};
%                where f_in is a single input value, and y_out is the
%                computed output value.
%             Example: To compute the step response:
%                   for j=1:100
%                      y(j) = iirdf1(1);
%                   end
%
% See also: iirdf2
%
% Version:  1.0
% Author:   D. Rowell   11/18/07
% -------------------------------------------------------------------------
% 
function y_n = iirdf1(f_n,B,A)
persistent f_register y_register Bx Ax N
%
% The following is initialization, and is executed once
%
if (ischar(f_n) && strcmp(f_n,'initial'))
    N = length(A);
    Ax = A;
    Bx = B;
    f_register = zeros(1,N);
    y_register = zeros(1,N);
    y_n = 0;
    else
% Filtering:  (Note that a Direct Form I filter needs two shift registers.)
    x = 0; y = 0;
    for J = N:-1:2
        y_register(J) = y_register(J-1); % Move along the shift register
        f_register(J) = f_register(J-1); 
        y = y - Ax(J)*y_register(J);
        x = x + Bx(J)*f_register(J);
    end
    x = x + Bx(1)*f_n;
    y_n = y + x;
    f_register(1) = f_n;
    y_register(1) = y_n;
end
