diff options
Diffstat (limited to 'macros/medfilt1.sci')
| -rw-r--r-- | macros/medfilt1.sci | 345 |
1 files changed, 345 insertions, 0 deletions
diff --git a/macros/medfilt1.sci b/macros/medfilt1.sci new file mode 100644 index 0000000..b81a03a --- /dev/null +++ b/macros/medfilt1.sci @@ -0,0 +1,345 @@ +function y = medfilt1(x, varargin) + // 1D median filtering + // + // Calling sequence + // y = medfilt1(x) + // y = medfilt1(x, n) + // y = medfilt1(x, n, dim) + // y = medfitl1(__, nanflag, padding) + // + // Description + // y = medfilt1(x) + // Applies a 3rd order 1-dimensional median filter to input x along the + // first non-zero dimension. The function appropriately pads the signal + // with zeros at the endings. For a segment, a median is calculated as + // the middle value (average of two middle values) for odd number + // number (even number) of data points. + // y = medfilt1(x,n) + // Applies a nth order 1-dimensional median filter. + // y = medfilt1(x,n,dim) + // Applies the median filter along the n-th dimension + // y = medfilt1(__, nanflag, padding) + // nanflag specifies how NaN values are treated. padding specifies the + // type of filtering to be performed at the signal edges. + // + // Parameters + // x: int | double + // Input signal. + // n: positive integer scalar + // Filter order. + // Defaults to 3.The order of the median filter. Must be less than + // (length of the signal) where signals are 1D vectors along the + // dimension of x to be filtered + // dim: positive integer scalar + // Dimension to filter along. + // Defaults to first non-singleton dimension of x + // nanflag: 'includenan' (default) | 'omitnan' + // NaN condition. + // * includenan: Filtering such that the median of any segment + // containing a NaN is also a NaN. + // * omitnan: Filtering with NaNs omitted in each segment. If a segment + // contains all NaNs, the result is NaN + // y: int | double + // The filtered signal. + // y has the same size as x + // + // Examples + // 1) Noise supression using median filtering + // fs = 1e3; + // t = 1:1/fs:1; + // s = sin(2*%pi*2*t)+ cos(2*%pi*5*t); + // // Adding noise + // x = s + 0.1*randn(size(s)); + // y = medfilt1(x); + // + // See also + // filter | hampel | median | sgolayfilt + // + // Authors + // Ayush Baid + + + + // ************************************************************************* + // Checking number of arguments + // ************************************************************************* + [numOutArgs, numInArgs] = argn(0); + + if numInArgs<1 | numInArgs>5 then + msg = "medfilt1: Wrong number of input argument; 1-5 expected"; + error(77, msg); + end + if numOutArgs~=1 then + msg = "medfilt1: Wrong number of output argument; 1 expected"; + error(78, msg); + end + + + + // ************************************************************************* + // Parsing input arguments + // ************************************************************************* + + // * Parsing x * + temp = x(:); + if type(temp)~=1 & type(temp)~=8 then + msg = "medfilt1: Wrong type for argument #1 (x): Int/double expected" + error(53, msg); + end + + + // * Parsing nanflag and padding * + // Getting all the string arguments + stringIndices = list(); + for i=1:length(varargin); + e = varargin(i); + if type(e)==10 then + stringIndices($+1)=i; + end + end + + nanflag = %f; // 0->includenan (default); 1->omitnan + padflag = %t; // 1->zeropad (default); 0->truncate + if ~isempty(stringIndices) then + // checking for 'omitnan' + if or(strcmpi(varargin(stringIndices), 'omitnan')) then + nanflag = %t; + end + + // checking for 'truncate' + if or(strcmpi(varargin(stringIndices), 'truncate')) then + padflag = %f; + end + varargin(stringIndices) = []; + end + + + // setting default value for n and dim + n = 3; + dim = 1; + L = length(size(x)); + for i=1:L + if size(x, i)>1 then + dim = i; + end + end + + // * Parsing n and dim * + if length(varargin)==1 then + if ~isempty(varargin(1)) then + n = varargin(1); + end + elseif length(varargin)==2 then + if ~isempty(varargin(1)) then + n = varargin(1); + end + if ~isempty(varargin(2)) then + dim = varargin(2); + end + else + msg = "medfilt1: Wrong type of input arguments; Atmost 3 numerical input expected"; + error(53, msg); + end + + // check on n + if length(n)~=1 then + msg = "medfilt1: Wrong size for argument #2 (n): Scalar expected"; + error(60,msg); + end + + if type(n)~=1 & type(n)~=8 then + msg = "medfilt1: Wrong type for argument #2 (n): Natural number expected"; + error(53,msg); + end + + if n~=round(n) | n<=0 then + msg = "medfilt1: Wrong type for argument #2 (n): Natural number expected"; + error(53,msg); + end + + if ~isreal(n) then + msg = "medfilt1: Wrong type for argument #2 (n): Real scalar expected"; + error(53,msg); + end + + // check on dim + if length(dim)~=1 then + msg = "medfilt1: Wrong size for argument #3 (dim): Scalar expected"; + error(60,msg); + end + + if type(dim)~=1 & type(dim)~=8 then + msg = "medfilt1: Wrong type for argument #3 (dim): Natural number expected"; + error(53,msg); + end + + if dim~=round(dim) | dim<=0 then + msg = "medfilt1: Wrong type for argument #3 (dim): Natural number expected"; + error(53,msg); + end + + if ~isreal(dim) then + msg = "medfilt1: Wrong type for argument #3 (dim): Real scalar expected"; + error(53,msg); + end + + + // ************************************************************************* + // Processing for median filtering column by column + // ************************************************************************* + + inp_size = size(x); + + + // Permuting x to bring the dimension to be acted upon as the first dimesnion + perm_vec = [2:dim, 1, dim+1:length(inp_size)]; + reverse_perm_vec = [dim, 1:dim-1, dim+1:length(inp_size)]; + x = permute(x, perm_vec); + + size_vec = size(x); + + y = x; // just initialization + + for i=1:prod(size_vec(2:$)) + temp = medfilt_colvector(x(:,i), n, padflag, nanflag); + y(:,i) = temp; + end + + + + y = permute(y, reverse_perm_vec); + + +endfunction + +function med = medfilt_colvector(x, n, zeropadflag, nanflag) + // Performs median filtering (of order n) on a column vector (x) + // zeropadflag -> zero pad instead of truncation + // nanflag -> discard all blocks containing nan, else do not consider nan values + + med = zeros(size(x,1),1); + disp('here1'); + + + // ** zero pad the signal ** + pad_length = floor(n/2); // padding on a size + x = [zeros(pad_length,1); x; zeros(pad_length,1)]; + + nx = length(x); + + // Arrange data in blocks + top_row = 1:(nx-n); + + idx = zeros(n,length(top_row)); + + for i=1:n + idx(i,:) = top_row + (i-1); + end + + blocks = matrix(x(idx), size(idx)); + + + if nanflag then + disp('here2'); + med = median(blocks, 1)'; + + // set result of all the blocks containing nan to nan + nanpresent = or(isnan(blocks), 1); + med(nanpresent) = %nan; + else + disp('here3'); + // we have to neglect nans + sorted_blocks = gsort(blocks, 'r', 'i'); + + // get the count of non-nan elements + num_elems = n - sum(isnan(sorted_blocks), 1); + + // find the median + offset = (0:size(blocks,2)-1)*size(blocks,1); + idx1 = offset+ceil(num_elems/2); + idx2 = offset+ceil((num_elems/2)+0.25); + + + // temporarily setting idx1 to 1 so as to not give errors in median calc. + // Will later replace values at such indices with Nan + idx1(idx1==0)=1; + med = (sorted_blocks(idx1) + sorted_blocks(idx2))./2; + + med(idx1==0) = %nan; + end + + if ~zeropadflag then + // ** recalculate boundary blocks with truncation truncate at the boundaries ** + + // divide the input signal into 3 parts; 1st and last part have truncation + for i=ceil(n/2):n + // ** first part ** + block = x(1:i); + + // * median calc for a block * + if nanflag then + med(i-ceil(n/2)+1) = median(block, 1); + + // set result of all the blocks containing nan to nan + nanpresent = or(isnan(block), 1); + if nanpresent then + med(i-ceil(n/2)+1) = %nan; + end + else + // we have to neglect nans + sorted_block = gsort(block, 'r', 'i'); + + // get the count of non-nan elements + num_elems = length(block) - sum(isnan(sorted_block), 1); + + // find the median + idx1 = ceil(num_elems/2); + idx2 = ceil(num_elems/2+0.25); + + + // temporarily setting idx1 to 1 so as to not give errors in median calc. + // Will later replace values at such indices with Nan + if idx1==0 then + med(i-ceil(n/2)+1) = %nan; + else + med(i-ceil(n/2)+1) = (sorted_block(idx1, :)+sorted_block(idx2, :))./2; + end + end + + + // ** last part ** + block = x($:-1:$-i); + + // * median calc for a block * + if nanflag then + med($+ceil(n/2)-i) = median(block, 1); + + // set result of all the blocks containing nan to nan + nanpresent = or(isnan(block), 1); + if nanpresent then + med($-ceil(n/2)+i) = %nan; + end + med($+ceil(n/2)-i) = %nan; + else + // we have to neglect nans + sorted_block = gsort(block, 'r', 'i'); + + // get the count of non-nan elements + num_elems = length(block) - sum(isnan(sorted_block), 1); + + // find the median + idx1 = ceil(num_elems/2); + idx2 = ceil(num_elems/2+0.25); + + // temporarily setting idx1 to 1 so as to not give errors in median calc. + // Will later replace values at such indices with Nan + if idx1==0 then + med($+ceil(n/2)-i) = %nan; + else + med($+ceil(n/2)-i) = (sorted_block(idx1) + sorted_block(idx2))./2; + end + end + end + end + +endfunction |