function [y]=MMP_recur_filt(x,tau);
%function [y]=MMP_recur_filt(x,tau)
%routine to apply Mudge & Lueck's version of the 
%  single pole filter recursive filter  JAOTech, 11, 825-836, 1994
% x are the input data
% tau is lag in scans;
% y are the filtered output data

%  Note discrepancy in Mudge & Lueck's notation and Matlab's filter.m notation
%  for a and b.  (They are reversed.)

%Mudge & Lueck's expression for the bilinear transform recursive filter is
%
%  y(n)=A*(x(n)+x(n-1)) + B*y(n-1)
%
% A = 1/(1+(2*tau/Ts)); B = 1-2A;
% 

%from Matlab's filter.m help, y=filter(b,a,x) outputs the following:
%   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)
% 

% with tau expressed in scans, the coefficients become:
%   note a = A,  b = B here to be consistent with Matlab

b(1)=(1/(1+2*tau));
b(2)=b(1);
a(1)=1;
a(2)=-(1-2*b(1));

y=filter(b,a,x-x(1));

y=y+x(1);

return
