Home > chronux_2_00 > spikesort > utility > datatools > filterm.m

filterm

PURPOSE ^

FILTERM One-dimensional digital filter (memory efficient).

SYNOPSIS ^

function [Y, carry] = filterm(B, A, X, carry)

DESCRIPTION ^

FILTERM           One-dimensional digital filter (memory efficient).
   Y = FILTERM(B,A,X) filters the data in vector X with the filter
   described by vectors A and B to create the filtered data Y.  The
   filter is a "Direct Form II Transposed" implementation of the standard
   difference equation: 

   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)
  
   If a(1) is not equal to 1, FILTER normalizes the filter
   coefficients by a(1).  

   If X is a matrix, FILTERM operates along the columns.  X is not
   allowed to be of dimension greater than 2.

   [Y,Zf] = FILTERM(B,A,X,Zi) gives access to initial and final
   conditions, Zi and Zf, of the delays.  Zi is a vector of length Q if
   X is a vector or of length Q x N if X is an M x N matrix, where
   Q = MAX(LENGTH(A),LENGTH(B))-1.  The returned Zf is always of type
   double, although Zi can be of any numeric data type.

   FILTERM is similar to Matlab's built-in FILTER except that it operates
   on non-double data, avoiding the memory penalty of converting an
   entire data array to double precision.  (The only difference in
   calling the two functions is that FILTERM will not operate on arrays
   of dimension greater than 2).

   If X is of type DOUBLE, INT8, INT16 or INT32, Y will be of the same
   type, with values clipped to the range of the original data type.  If
   X is of an unsigned integer type (UINT8, UINT16, UINT32), Y will be of
   the corresponding signed type.  For example, if X is of type UINT8, Y
   will be of type INT8 and values larger than [-128,127] will be
   clipped.  If the data X uses the full range of an unsigned data type,
   it should be first converted to a data type with larger dynamic range.

   Note that if A = 1, it is more efficient to use FILTERZM.

   See also FILTER, FILTERZM.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Y, carry] = filterm(B, A, X, carry)
0002 %FILTERM           One-dimensional digital filter (memory efficient).
0003 %   Y = FILTERM(B,A,X) filters the data in vector X with the filter
0004 %   described by vectors A and B to create the filtered data Y.  The
0005 %   filter is a "Direct Form II Transposed" implementation of the standard
0006 %   difference equation:
0007 %
0008 %   a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
0009 %                           - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
0010 %
0011 %   If a(1) is not equal to 1, FILTER normalizes the filter
0012 %   coefficients by a(1).
0013 %
0014 %   If X is a matrix, FILTERM operates along the columns.  X is not
0015 %   allowed to be of dimension greater than 2.
0016 %
0017 %   [Y,Zf] = FILTERM(B,A,X,Zi) gives access to initial and final
0018 %   conditions, Zi and Zf, of the delays.  Zi is a vector of length Q if
0019 %   X is a vector or of length Q x N if X is an M x N matrix, where
0020 %   Q = MAX(LENGTH(A),LENGTH(B))-1.  The returned Zf is always of type
0021 %   double, although Zi can be of any numeric data type.
0022 %
0023 %   FILTERM is similar to Matlab's built-in FILTER except that it operates
0024 %   on non-double data, avoiding the memory penalty of converting an
0025 %   entire data array to double precision.  (The only difference in
0026 %   calling the two functions is that FILTERM will not operate on arrays
0027 %   of dimension greater than 2).
0028 %
0029 %   If X is of type DOUBLE, INT8, INT16 or INT32, Y will be of the same
0030 %   type, with values clipped to the range of the original data type.  If
0031 %   X is of an unsigned integer type (UINT8, UINT16, UINT32), Y will be of
0032 %   the corresponding signed type.  For example, if X is of type UINT8, Y
0033 %   will be of type INT8 and values larger than [-128,127] will be
0034 %   clipped.  If the data X uses the full range of an unsigned data type,
0035 %   it should be first converted to a data type with larger dynamic range.
0036 %
0037 %   Note that if A = 1, it is more efficient to use FILTERZM.
0038 %
0039 %   See also FILTER, FILTERZM.
0040 
0041 % Based on from TMW's FILTER.
0042 %  The following commented code tests FILTERM against FILTER:
0043 %   (also, monitor memory usage to see memory advantage)
0044 % X = int8(randn(2^22,2)*32);   [B,A] = butter(6,0.1);
0045 % disp('Running FILTER ...');   drawnow;  tic;
0046 % Y  = int8(filter(B, A, double(X)));  t(1) = toc;
0047 % disp('Running FILTERM ...');  drawnow;  tic;
0048 % YM = filterm(B, A, X);               t(2) = toc;
0049 % check = ceil(rand(1e6,1)*length(X));
0050 % mse = sqrt(sum(double(Y(check))-double(YM(check)).^2));
0051 % printf('FILTER: %4.3f sec   FILTERM: %4.3f sec     (APPROX) MSE: %5.3f', t(1), t(2), mse);
0052 
0053 
0054 % # elements to process at once: for some reason ~2^11-2^13 seems to work
0055 % best (probably cpu cache, although the range is similar on computers
0056 % with different cache sizes ...) & exact powers of 2 sometimes cause
0057 % hiccups ... so heuristically, 2000 seems reasonable
0058 chunksize = 2000; 
0059 
0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Check Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0061 % filter coefficients check
0062 if (A(1) ~= 1),  A = A ./ A(1);  end;   % normalize b/c FILTER does
0063 
0064 % input data check
0065 if (ndims(X) > 2), error('FILTERM is undefined for arrays of dimension greater than 2.');  end;
0066 if (isvectord(X) == 2), T = 1;  X = X';  else  T = 0;  end;  % col vect for now
0067 [M,N] = size(X);
0068 
0069 % initial conditions check
0070 Q = max(length(A), length(B)) - 1;
0071 if (nargin > 3)
0072     if (isvectord(carry)), carry = carry(:);  end;  % force column vector
0073     [Qc,Nc] = size(carry);
0074     if (Q~=Qc || N~=Nc)
0075         error('Initial conditions array is not of the correct size.');
0076     end
0077     carry = double(carry);
0078 else  carry = zeros(Q,N);
0079 end
0080 
0081 % data type check
0082 switch (class(X)),
0083     case 'double',              convert = @double;
0084     case {'uint8','int8'},    convert = @int8;
0085     case {'uint16','int16'},  convert = @int16;
0086     case {'uint32','int32'},  convert = @int32;
0087     otherwise, error('X must be a numeric data type.');
0088 end
0089 
0090 
0091 %%%%%%%%%%%%%%%%%%%%%%%%%% Do the Filtering %%%%%%%%%%%%%%%%%%%%%%%%%%
0092 if (isa(X, 'double') || M < chunksize),  % no point chunking if double
0093     [Y,carry] = filter(B, A, double(X), carry);
0094     Y = feval(convert, Y);
0095 else
0096     % Perform the filtering in chunks so that only a portion of the data
0097     % needs to be expanded to double at any given time.
0098     Y = repmat(feval(convert,0), size(X));
0099     for start = 1:chunksize:M
0100         finish = min(M, start+chunksize-1);
0101         [chunk,carry] = filter(B,A,double(X(start:finish,:)),carry);
0102         Y(start:finish,:) = feval(convert, chunk);
0103     end
0104 end
0105 
0106 
0107 %%%%%%%%%%%%%%%%%%%%%%%%%% Clean Up Outputs %%%%%%%%%%%%%%%%%%%%%%%%%%
0108 if (T),  Y = Y';  end;   % if input was row vect, keep same orientation
0109 if (nargout < 2), clear carry; end;

Generated on Fri 15-Aug-2008 11:35:42 by m2html © 2003