


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.

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;