


FILTERZM Zero-phase 1-D digital filter (memory efficient).
Y = FILTERZ(B, X) filters the data in vector X with the FIR filter
given by the vector B to create the filtered data Y, shifting the
filter to give zero-phase delay. The process is equivalent to
convolution with the time-reversed and shifted filter B:
For P = ceil((length(B)+1)/2),
y(n) = ... + x(n-1)*b(P+1) + x(n)*b(P) + x(n+1)*b(P-1) + ...
If X is a matrix, FILTERZM operates along the columns. X is not
allowed to be of dimension greater than 2.
FILTERZM(B, X) is similar to FILTER2(B,X,'same'), except that FILTERZM
operates on non-double data, avoiding the memory penalty of converting
an entire data array to double precision. It is analogous in this
respect to FILTERM and follows the same datatype conventions as that
function. However, by restricting usage only to FIR filters, FILTERZM
provides zero-phase delay filtering and is faster than FILTERM(B,1,X).
See also FILTERM, FILTER2.

0001 function Y = filterzm(B, X) 0002 %FILTERZM Zero-phase 1-D digital filter (memory efficient). 0003 % Y = FILTERZ(B, X) filters the data in vector X with the FIR filter 0004 % given by the vector B to create the filtered data Y, shifting the 0005 % filter to give zero-phase delay. The process is equivalent to 0006 % convolution with the time-reversed and shifted filter B: 0007 % For P = ceil((length(B)+1)/2), 0008 % y(n) = ... + x(n-1)*b(P+1) + x(n)*b(P) + x(n+1)*b(P-1) + ... 0009 % 0010 % If X is a matrix, FILTERZM operates along the columns. X is not 0011 % allowed to be of dimension greater than 2. 0012 % 0013 % FILTERZM(B, X) is similar to FILTER2(B,X,'same'), except that FILTERZM 0014 % operates on non-double data, avoiding the memory penalty of converting 0015 % an entire data array to double precision. It is analogous in this 0016 % respect to FILTERM and follows the same datatype conventions as that 0017 % function. However, by restricting usage only to FIR filters, FILTERZM 0018 % provides zero-phase delay filtering and is faster than FILTERM(B,1,X). 0019 % 0020 % See also FILTERM, FILTER2. 0021 0022 % Based on from TMW's FILTER. 0023 % The following commented code tests FILTERMZ against FILTER: 0024 % (also, monitor memory usage to see memory advantage) 0025 % X = int8(rand(2^22,2)*64); B = fir1(6, 0.1); 0026 % disp('Running FILTER ...'); drawnow; tic; 0027 % Y = int8(filter(B, 1, double(X))); t(1) = toc; 0028 % disp('Running FILTERZM ...'); drawnow; tic; 0029 % YM = filterzm(B, X); t(2) = toc; 0030 % check = ceil(rand(1e6,1)*length(X)); 0031 % mse = sqrt(sum(double(Y(check))-double(YM(check)).^2)); 0032 % printf('FILTER: %4.3f sec FILTERZM: %4.3f sec (APPROX) MSE: %5.3f', t(1), t(2), mse); 0033 0034 0035 % # elements to process at once: for some reason ~2^11-2^13 seems to work 0036 % best (probably cpu cache, although the range is similar on computers 0037 % with different cache sizes ...) & exact powers of 2 sometimes cause 0038 % hiccups ... so heuristically, 2000 seems reasonable 0039 chunksize = 2000; 0040 0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Check Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 % input data check 0043 if (ndims(X) > 2), error('FILTERZM is undefined for arrays of dimension greater than 2.'); end; 0044 if (isvectord(X) == 2), T = 1; X = X'; else T = 0; end; % col vect for now 0045 [M,N] = size(X); 0046 0047 % Massage filter into proper shape 0048 if (~isvectord(B)), error('The FIR filter B must be 1-dimensional.'); end; 0049 if (isvectord(B) == 2), B = B'; end; % force column vector 0050 B = flipud(B); 0051 P = (length(B)-1)/2; Pt = floor(P); Pb = ceil(P); % top/bot filter overhang 0052 0053 % data type check 0054 switch (class(X)), 0055 case 'double', convert = @double; 0056 case {'uint8','int8'}, convert = @int8; 0057 case {'uint16','int16'}, convert = @int16; 0058 case {'uint32','int32'}, convert = @int32; 0059 otherwise, error('X must be a numeric data type.'); 0060 end 0061 0062 0063 %%%%%%%%%%%%%%%%%%%%%%%%%% Do the Filtering %%%%%%%%%%%%%%%%%%%%%%%%%% 0064 if (isa(X, 'double') || M < chunksize), % no point chunking if double 0065 Y = conv2(double(X), B, 'same'); 0066 Y = feval(convert, Y); 0067 else 0068 % Perform the filtering in chunks so that only a portion of the data 0069 % needs to be expanded to double at any given time. 0070 Y = repmat(feval(convert,0), size(X)); 0071 X = [repmat(feval(convert,0), [Pt+1, N]); X; repmat(feval(convert,0), [Pb+1, N])]; 0072 for start = 1:chunksize:M 0073 finish = min(M, start+chunksize-1); 0074 chunk = conv2(double((X(start:(finish+Pt+Pb+2),:))),B,'same'); 0075 Y(start:finish,:) = feval(convert, chunk((Pt+2):(end-Pb-1),:)); 0076 end 0077 end 0078 0079 0080 %%%%%%%%%%%%%%%%%%%%%%%%%% Clean Up Outputs %%%%%%%%%%%%%%%%%%%%%%%%%% 0081 if (T), Y = Y'; end; % if input was row vect, keep same orientation