Description of Smooth

0001 function smoothed = Smooth(data,smooth,varargin)
0002 
0003 %Smooth - Smooth using a Gaussian kernel.
0004 %
0005 %  USAGE
0006 %
0007 %    smoothed = Smooth(data,smooth,<options>)
0008 %
0009 %    data           data to smooth
0010 %    smooth         vertical and horizontal kernel sizes:
0011 %                    - gaussian: standard deviations [Sv Sh] (optionally,
0012 %                      kernel size W can also be provided as [Sv Sh Wv Wh])
0013 %                    - rectangular/triangular: window half size [Nv Nh]
0014 %                   (in number of samples, 0 = no smoothing)
0015 %    <options>      optional list of property-value pairs (see table below))
0016 %
0017 %    =========================================================================
0018 %     Properties    Values
0019 %    -------------------------------------------------------------------------
0020 %     'type'        two letters (one for X and one for Y) indicating which
0021 %                   coordinates are linear ('l') and which are circular ('c')
0022 %                   - for 1D data, only one letter is used (default 'll')
0023 %     'kernel'      either 'gaussian' (default), 'rectangular' (running
0024 %                   average), or 'triangular' (weighted running average)
0025 %    =========================================================================
0026 %
0027 
0028 % Copyright (C) 2004-2016 by Michaël Zugaro, 2013 Nicolas Maingret
0029 %
0030 % This program is free software; you can redistribute it and/or modify
0031 % it under the terms of the GNU General Public License as published by
0032 % the Free Software Foundation; either version 3 of the License, or
0033 % (at your option) any later version.
0034 
0035 maxSize = 10001;
0036 
0037 if nargin < 2,
0038     error('Incorrect number of parameters (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).');
0039 end
0040 
0041 vector = isvector(data);
0042 matrix = (~vector & length(size(data)) == 2);
0043 if ~vector & ~matrix,
0044     error('Smoothing applies only to vectors or matrices (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).');
0045 end
0046 
0047 % Vectors must be 'vertical'
0048 if size(data,1) == 1,
0049     data = data';
0050 end
0051 
0052 % Default values
0053 kernel = 'gaussian';
0054 if vector, type = 'l'; else type = 'll'; end
0055 
0056 % If Sh = Sv = 0, no smoothing required
0057 if all(smooth==0),
0058     smoothed = data;
0059     return
0060 end
0061 
0062 % Parse parameter list
0063 for i = 1:2:length(varargin),
0064     if ~ischar(varargin{i}),
0065         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).']);
0066     end
0067     switch(lower(varargin{i})),
0068         case 'type',
0069             type = lower(varargin{i+1});
0070             if (vector && ~isastring(type,'c','l')) || (~vector && ~isastring(type,'cc','cl','lc','ll')),
0071                 error('Incorrect value for property ''type'' (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).');
0072             end
0073         case 'kernel',
0074             kernel = lower(varargin{i+1});
0075             if ~isastring(kernel,'gaussian','rectangular','triangular'),
0076                 error('Incorrect value for property ''kernel'' (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).');
0077             end
0078         otherwise,
0079             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).']);
0080 
0081   end
0082 end
0083 
0084 % Check kernel parameters
0085 if matrix && length(smooth) == 1,
0086     % For 2D data, providing only one value S for the std is interpreted as Sh = Sv = S
0087     smooth = [smooth smooth];
0088 end
0089 if strcmp(kernel,'gaussian'),
0090     if ~isdvector(smooth,'>=0') | (vector & ~ismember(length(smooth),[1 2])) | (matrix & ~ismember(length(smooth),[2 4])),
0091         error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).');
0092     end
0093 else
0094     if ~isdvector(smooth,'>=0') | (vector & length(smooth) ~= 1) | (matrix & length(smooth) ~= 2),
0095         error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help Smooth">Smooth</a>'' for details).');
0096     end
0097 end
0098 
0099 [vSize,hSize] = size(data);
0100 
0101 % Build kernels
0102 if strcmp(kernel,'rectangular'),
0103 
0104     % Rectangular kernel (running average)
0105     % 1) Vertical kernel
0106     vKernelSize = 2*smooth(1)+1;
0107     vKernel = ones(vKernelSize,1);
0108     vKernel = vKernel/sum(vKernel);
0109     % 2) Horizontal kernel
0110     if ~vector,
0111         hKernelSize = 2*smooth(2)+1;
0112         hKernel = ones(hKernelSize,1);
0113         hKernel = hKernel/sum(hKernel);
0114     end
0115 
0116 elseif strcmp(kernel,'triangular'),
0117 
0118     % Triangular kernel
0119     % 1) Vertical kernel
0120     vKernelSize = 2*smooth(1)+1;
0121     vKernel = triang(vKernelSize);
0122     vKernel = vKernel/sum(vKernel);
0123     % 2) Horizontal kernel
0124     if ~vector,
0125         hKernelSize = 2*smooth(2)+1;
0126         hKernel = ones(hKernelSize,1);
0127         hKernel = hKernel/sum(hKernel);
0128     end
0129 
0130 else
0131 
0132     % Gaussian kernel
0133     % 1) Vertical kernel
0134     if vector && length(smooth) == 2,
0135         vKernelSize = smooth(2);
0136     elseif matrix && length(smooth) == 4,
0137         vKernelSize = smooth(3);
0138     else
0139         if vSize > maxSize, warning(['Default kernel too large; using ' int2str(maxSize) ' points.']); end
0140         vKernelSize = min([vSize maxSize]);
0141     end
0142     r = (-vKernelSize:vKernelSize)'/vKernelSize;
0143     vKernelStdev = smooth(1)/vKernelSize;
0144     vKernel = exp(-r.^2/(vKernelStdev+eps)^2/2);
0145     vKernel = vKernel/sum(vKernel);
0146     % 2) Horizontal kernel
0147     if ~vector,
0148         if length(smooth) == 4,
0149             hKernelSize = smooth(4);
0150         else
0151             if hSize > maxSize, warning(['Default kernel too large; using ' int2str(maxSize) ' points.']); end
0152             hKernelSize = min([hSize maxSize]);
0153         end
0154         r = (-hKernelSize:hKernelSize)/hKernelSize;
0155         hKernelStdev = smooth(2)/hKernelSize;
0156         hKernel = exp(-r.^2/(hKernelStdev+eps)^2/2);
0157         hKernel = hKernel/sum(hKernel);
0158     end
0159     
0160 end
0161 
0162 if vector,
0163     % Vector smoothing
0164     % Prepend/append data to limit edge effects
0165     if strcmp(type,'l'),
0166         % For linear data, flip edge data
0167         % top = 2*data(1)-flipud(data(1:vKernelSize));
0168         % bottom = 2*data(end)-flipud(data(end-vKernelSize+1:end));
0169         top = flipud(data(1:vKernelSize));
0170         bottom = flipud(data(end-vKernelSize+1:end));
0171     else
0172         % For circular data, wrap edge data
0173         top = data(end-vKernelSize+1:end);
0174         bottom = data(1:vKernelSize);
0175     end
0176     data = [top;data;bottom];
0177     % Convolve (and return central part)
0178     tmp = conv(vKernel,data);
0179     n = size(tmp,1);
0180     d = n - vSize;
0181     start = d/2+1;
0182     stop = start + vSize - 1;
0183     smoothed = tmp(start:stop,:);
0184 else
0185     % Matrix smoothing
0186     % Convolve
0187     if smooth(1) == 0,
0188         % Smooth only across columns (Sv = 0)
0189         % Prepend/append data to limit edge effects
0190         if strcmp(type(1),'l'),
0191             % For linear data, flip edge data
0192             % left = 2*repmat(data(:,1),1,hKernelSize)-fliplr(data(:,1:hKernelSize));
0193             % right = 2*repmat(data(:,end),1,hKernelSize)-fliplr(data(:,end-hKernelSize+1:end));
0194             left = fliplr(data(:,1:hKernelSize));
0195             right = fliplr(data(:,end-hKernelSize+1:end));
0196         else
0197             % For circular data, wrap edge data
0198             left = data(:,end-hKernelSize+1:end);
0199             right = data(:,1:hKernelSize);
0200         end
0201         data = [left data right];
0202         for i = 1:size(data,1),
0203             tmp = conv(hKernel,data(i,:));
0204             n = size(tmp,2);
0205             d = n - hSize;
0206             start = d/2+1;
0207             stop = start + hSize - 1;
0208             smoothed(i,:) = tmp(:,start:stop);
0209         end
0210     elseif smooth(2) == 0,
0211         % Smooth only across lines (Sh = 0)
0212         % Prepend/append data to limit edge effects
0213         if strcmp(type(2),'l'),
0214             % For linear data, flip edge data
0215             % top = 2*repmat(2*data(1,:),vKernelSize,1)-flipud(data(1:vKernelSize,:));
0216             % bottom = 2*repmat(2*data(end,:),vKernelSize,1)-flipud(data(end-vKernelSize+1:end,:));
0217             top = flipud(data(1:vKernelSize,:));
0218             bottom = flipud(data(end-vKernelSize+1:end,:));
0219         else
0220             % For circular data, wrap edge data
0221             bottom = data(1:vKernelSize,:);
0222             top = data(end-vKernelSize+1:end,:);
0223         end
0224         data = [top;data;bottom];
0225         for i = 1:size(data,2),
0226             tmp = conv(vKernel,data(:,i));
0227             n = size(tmp,1);
0228             d = n - vSize;
0229             start = d/2+1;
0230             stop = start + vSize - 1;
0231             smoothed(:,i) = tmp(start:stop);
0232         end
0233     else
0234         % Smooth in 2D
0235         % Prepend/append data to limit edge effects
0236         if strcmp(type(2),'l'),
0237             % For linear data, flip edge data
0238               % top = 2*repmat(2*data(1,:),vKernelSize,1)-flipud(data(1:vKernelSize,:));
0239               % bottom = 2*repmat(2*data(end,:),vKernelSize,1)-flipud(data(end-vKernelSize+1:end,:));
0240             top = flipud(data(1:vKernelSize,:));
0241             bottom = flipud(data(end-vKernelSize+1:end,:));
0242         else
0243             % For circular data, wrap edge data
0244             bottom = data(1:vKernelSize,:);
0245             top = data(end-vKernelSize+1:end,:);
0246         end
0247         data = [top;data;bottom];
0248         if strcmp(type(1),'l'),
0249             % For linear data, flip edge data
0250             % left = 2*repmat(2*data(:,1),1,hKernelSize)-fliplr(data(:,1:hKernelSize));
0251               % right = 2*repmat(2*data(:,end),1,hKernelSize)-fliplr(data(:,end-hKernelSize+1:end));
0252             left = fliplr(data(:,1:hKernelSize));
0253             right = fliplr(data(:,end-hKernelSize+1:end));
0254         else
0255             % For circular data, wrap edge data
0256             left = data(:,end-hKernelSize+1:end);
0257             right = data(:,1:hKernelSize);
0258         end
0259         data = [left data right];
0260         tmp = conv2(vKernel,hKernel,data,'same');
0261         n = size(tmp,1);
0262         d = n - vSize;
0263         vStart = d/2+1;
0264         vStop = vStart + vSize - 1;
0265         n = size(tmp,2);
0266         d = n - hSize;
0267         hStart = d/2+1;
0268         hStop = hStart + hSize - 1;
0269         smoothed = tmp(vStart:vStop,hStart:hStop);
0270     end
0271 end
Smooth

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
