medfilt

 medfilt -- Median-filtering.
  medfilt(theData, theFilterLength, theChunkLength) applies
   a centered median-filter of theFilterLength (default = 5)
   to the columns of theData.  The computations are segmented
   by the ChunkLength (default is the full height of the data).
   Zero-slope end-conditions on the original data are assumed.
   (Note: a row-vector will be filtered as a row.)

0001 function theResult = medfilt(theData, theFilterLength, theChunkLength)
0002 % medfilt -- Median-filtering.
0003 %  medfilt(theData, theFilterLength, theChunkLength) applies
0004 %   a centered median-filter of theFilterLength (default = 5)
0005 %   to the columns of theData.  The computations are segmented
0006 %   by the ChunkLength (default is the full height of the data).
0007 %   Zero-slope end-conditions on the original data are assumed.
0008 %   (Note: a row-vector will be filtered as a row.)
0009  
0010 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0011 %  All Rights Reserved.
0012 %   Disclosure without explicit written consent from the
0013 %    copyright owner does not constitute publication.
0014  
0015 % Version of 14-May-1998 11:37:36.
0016 
0017 if nargin < 1, theData = 'demo'; end
0018 
0019 if isequal(theData, 'demo')
0020         help(mfilename)
0021         theData = fix(20 .* rand(8, 2))
0022         theFilterLength = 3
0023         theChunkLength = 5
0024         result = medfilt(theData, theFilterLength, theChunkLength)
0025         return
0026 end
0027 
0028 oldSize = size(theData);
0029 m = oldSize(1);
0030 n = prod(oldSize) ./ m;
0031 theData = reshape(theData, [m n]);
0032 if m == 1 & length(oldSize) == 2
0033         theData = reshape(theData, [n m]);
0034 end
0035 
0036 [m, n] = size(theData);
0037 
0038 if nargin < 2, theFilterLength = 5; end
0039 if nargin < 3, theChunkLength = m; end
0040 
0041 result = zeros(size(theData));
0042 
0043 if rem(theFilterLength, 2) == 0
0044         theFilterLength = theFilterLength + 1;
0045 end
0046 
0047 f = fix(theFilterLength./2);
0048 
0049 k = 0;
0050 while k < m
0051         kmin = k+1;
0052         kmax = min(k+theChunkLength, m);
0053         imin = max(1, kmin-f);
0054         imax = min(kmax+f, m);
0055         d = theData(imin:imax, :);
0056         r = mdfilt(d, theFilterLength);
0057         imax = imax - imin + 1;
0058         imin = 1;
0059         if kmax < m, imax = imax - f; end
0060         if kmin > 1, imin = imin + f; end
0061         result(kmin:kmax, :) = r(imin:imax, :);
0062         k = kmax;
0063 end
0064 
0065 result = reshape(result, oldSize);
0066 
0067 if nargout > 0
0068         theResult = result;
0069 else
0070         assignin('caller', 'ans', result)
0071         ans = result
0072 end
0073 
0074 % -------------------- mdfilt -------------------- %
0075 
0076 function theResult = mdfilt(theData, theFilterLength)
0077 
0078 % mdfilt -- Median filter.
0079 %  mdfilt(theData, theFilterLength) applies a centered
0080 %   median-filter of theFilterLength (odd-number) to the
0081 %   columns of theData, assuming zero-slope ends.  The
0082 %   default filter-length is 5.  The result is the same
0083 %   size as the data.  (Note: a row-vector is filtered
0084 %   as a row.)
0085 %  medfilt('demo') demonstrates itself.
0086  
0087 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0088 %  All Rights Reserved.
0089 %   Disclosure without explicit written consent from the
0090 %    copyright owner does not constitute publication.
0091  
0092 % Version of 11-May-1998 16:59:16.
0093 % Revised    12-May-1998 23:58:38.
0094 
0095 if nargin < 1, theData = 'demo'; end
0096 
0097 if isequal(theData, 'demo')
0098         help(mfilename)
0099         m = 6; n = 2;
0100         theData = fix(rand(m, n).*20)
0101         theFilterLength = 3
0102         theFilteredData = mdfilt(theData, theFilterLength)
0103         return
0104 end
0105 
0106 if nargin < 2, theFilterLength = 5; end
0107 
0108 % Use odd-length filter.
0109 
0110 if rem(theFilterLength, 2) == 0
0111         theFilterLength = theFilterLength + 1;
0112 end
0113 
0114 % Re-arrange for efficiency.
0115 
0116 oldSize = size(theData);
0117 m = oldSize(1); n = prod(oldSize)./m;
0118 if m == 1 & length(oldSize) == 2
0119         theData = theData.';
0120         [m, n] = size(theData);
0121 end
0122 theData = reshape(theData, [m n]);
0123 theData = theData.';   % Transpose to row-wise.
0124 [m, n] = size(theData);
0125 
0126 result = zeros(size(theData));
0127 
0128 % Create index-array.
0129 
0130 f = fix(theFilterLength./2);
0131 ind = [ones(1, f) (1:n) ones(1, f).*n];
0132 indices = zeros(theFilterLength, n);
0133 k = 1:n;
0134 for i = 1:theFilterLength
0135         indices(i, :) = ind(k);
0136         k = k + 1;
0137 end
0138 
0139 % Filter each row.
0140 
0141 x = zeros(theFilterLength, n);
0142 for i = 1:m
0143         x(:) = theData(i, indices);
0144         result(i, :) = median(x);
0145 end
0146 
0147 % Undo transpose and reshape.
0148 
0149 result = reshape(result.', oldSize);
0150 
0151 % Output.
0152 
0153 if nargout > 0
0154         theResult = result;
0155 else
0156         assignin('caller', 'ans', result)
0157         ans = result
0158 end