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

filterzm

PURPOSE ^

FILTERZM Zero-phase 1-D digital filter (memory efficient).

SYNOPSIS ^

function Y = filterzm(B, X)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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