nanmoving_average

PURPOSE ^

NANMOVING_AVERAGE Moving average ignoring NaN's.

SYNOPSIS ^

function [Y,Nsum] = nanmoving_average(X,F,DIM,INT)

DESCRIPTION ^

NANMOVING_AVERAGE   Moving average ignoring NaN's.

   Syntax:
     [Y,Nsum] = nanmoving_average(X,F,DIM,INT);

   Input:
     X   - Vector or matrix of finite elements.
     F   - Window semi-length. A positive scalar.
     DIM - If DIM=1: smooths the columns (default); elseif DIM=2 the rows.
     INT - If INT=0: do not interpolates NaN's (default); elseif INT=1 do
           interpolates.

   Output:
     Y    - Smoothed X elements with/without interpolated NaN's.
     Nsum - Number of not NaN's elements that fixed on the moving window.
            Provided to get a sum instead of a mean: Y.*Nsum.

   Description:
     Quickly smooths the vector X by averaging each element along with the
     2*F elements at its sides ignoring NaN's. The elements at the ends
     are also averaged but the extrems are left intact. With the windows
     size defined in this way, the filter has zero phase. If all the 2*F+1
     elements al NaN's, a NaN is return.

   Example:
      x = 2*pi*linspace(-1,1); 
      yn = cos(x) + 0.25 - 0.5*rand(size(x)); 
      yn([5 30 40:43]) = NaN;
      ys = nanmoving_average(yn,4);
      yi = nanmoving_average(yn,4,[],1);
      plot(x,yn,x,yi,'.',x,ys),  axis tight
      legend('noisy','smooth','without interpolation',4)

   See also FILTER, RECTWIN, NANMEAN and MOVING_AVERAGE, MOVING_AVERAGE2,
   NANMOVING_AVERAGE2 by Carlos Vargas

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Y,Nsum] = nanmoving_average(X,F,DIM,INT)
0002 %NANMOVING_AVERAGE   Moving average ignoring NaN's.
0003 %
0004 %   Syntax:
0005 %     [Y,Nsum] = nanmoving_average(X,F,DIM,INT);
0006 %
0007 %   Input:
0008 %     X   - Vector or matrix of finite elements.
0009 %     F   - Window semi-length. A positive scalar.
0010 %     DIM - If DIM=1: smooths the columns (default); elseif DIM=2 the rows.
0011 %     INT - If INT=0: do not interpolates NaN's (default); elseif INT=1 do
0012 %           interpolates.
0013 %
0014 %   Output:
0015 %     Y    - Smoothed X elements with/without interpolated NaN's.
0016 %     Nsum - Number of not NaN's elements that fixed on the moving window.
0017 %            Provided to get a sum instead of a mean: Y.*Nsum.
0018 %
0019 %   Description:
0020 %     Quickly smooths the vector X by averaging each element along with the
0021 %     2*F elements at its sides ignoring NaN's. The elements at the ends
0022 %     are also averaged but the extrems are left intact. With the windows
0023 %     size defined in this way, the filter has zero phase. If all the 2*F+1
0024 %     elements al NaN's, a NaN is return.
0025 %
0026 %   Example:
0027 %      x = 2*pi*linspace(-1,1);
0028 %      yn = cos(x) + 0.25 - 0.5*rand(size(x));
0029 %      yn([5 30 40:43]) = NaN;
0030 %      ys = nanmoving_average(yn,4);
0031 %      yi = nanmoving_average(yn,4,[],1);
0032 %      plot(x,yn,x,yi,'.',x,ys),  axis tight
0033 %      legend('noisy','smooth','without interpolation',4)
0034 %
0035 %   See also FILTER, RECTWIN, NANMEAN and MOVING_AVERAGE, MOVING_AVERAGE2,
0036 %   NANMOVING_AVERAGE2 by Carlos Vargas
0037 
0038 % Copyright 2006-2008  Carlos Vargas, nubeobscura@hotmail.com
0039 %    $Revision: 3.1 $  $Date: 2008/03/12 18:20:00 $
0040 
0041 %   Written by
0042 %   M. in S. Carlos Adrián Vargas Aguilera
0043 %   Physical Oceanography PhD candidate
0044 %   CICESE
0045 %   Mexico,  march 2008
0046 %
0047 %   nubeobscura@hotmail.com
0048 %
0049 %   Download from:
0050 %   http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
0051 %   tType=author&objectId=1093874
0052 
0053 %   2008 Mar. Use CUMSUM as RUNMEAN by Jos van der Geest, no more
0054 %   subfunctions.
0055 
0056 %% Error checking
0057 if ~nargin
0058   error('Nanmoving_average:Inputs','There are no inputs.')
0059 elseif nargin<2 || isempty(F)
0060  F = 0;
0061 end
0062 if F==0
0063  Y = X;
0064  return
0065 end
0066 F = round(F);
0067 ndim = ndims(X);
0068 if (ndim ~= 2)
0069  error('Nanmoving_average:Inputs','Input is not a vector or matrix.')
0070 end
0071 [N,M] = size(X);
0072 if nargin<3 || isempty(DIM)
0073  DIM = 1;
0074  if N == 1
0075   DIM = 2;
0076  end
0077 end
0078 if DIM == 2
0079  X = X.';
0080  [N,M] = size(X);
0081 end
0082 if 2*F+1>N
0083  warning('Nanmoving_average:Inputs',...
0084   'Window size must be less or equal as the number of elements.')
0085  Y = X;
0086  if DIM == 2
0087   Y = Y.';
0088  end
0089  return
0090 end
0091 if nargin<4 || isempty(INT)
0092  INT = 0;
0093 end
0094 
0095 %% Window width
0096 Wwidth = 2*F + 1;
0097 
0098 %% Smooth the edges but with the first and last element intact
0099 F2 = Wwidth - 2;
0100 Y1 =        X(     1:F2,:);
0101 Y2 = flipud(X(N-F2+1:N ,:));
0102 inan1 = isnan(Y1); 
0103 inan2 = isnan(Y2);
0104 Y1(inan1) = 0;
0105 Y2(inan2) = 0;
0106 Y1    = cumsum(Y1,1);    Y1    =    Y1(1:2:F2,:);
0107 Y2    = cumsum(Y2,1);    Y2    =    Y2(1:2:F2,:);
0108 inan1 = cumsum(inan1,1); inan1 = inan1(1:2:F2,:);
0109 inan2 = cumsum(inan2,1); inan2 = inan2(1:2:F2,:);
0110 Nsum1 = repmat((1:2:F2)',1,M);
0111 Nsum2 = repmat((1:2:F2)',1,M);
0112 Nsum1 = Nsum1 - inan1;
0113 Nsum2 = Nsum2 - inan2;
0114 Y1(Nsum1==0) = NaN;
0115 Y2(Nsum2==0) = NaN;
0116 Y1 = Y1./Nsum1;            
0117 Y2 = Y2./Nsum2;
0118 
0119 %% Recursive moving average method:
0120 nnan = ~isnan(X);
0121 X(~nnan) = 0;
0122 % Cumsum trick copied from RUNMEAN by Jos van der Geest  (12 mar 2008)
0123 Y = [zeros(F+1,M); X; zeros(F,M)];
0124 Y = cumsum(Y,1);
0125 Y = Y(Wwidth+1:end,:)-Y(1:end-Wwidth,:);
0126 Nsum = [zeros(F+1,M); nnan; zeros(F,M)];
0127 Nsum = cumsum(Nsum,1);
0128 Nsum = Nsum(Wwidth+1:end,:)-Nsum(1:end-Wwidth,:);
0129 Y = Y./Nsum;
0130 
0131 %% Sets the smoothed edges:
0132 Y(    1:F,:) =        Y1;
0133 Y(N-F+1:N,:) = flipud(Y2);
0134 Nsum(    1:F,:) = Nsum1;
0135 Nsum(N-F+1:N,:) = flipud(Nsum2);
0136 
0137 %% Do not interpolates:
0138 if ~INT
0139  Y(~nnan) = NaN;
0140 end
0141 
0142 %% Return correct size:
0143 if DIM == 2
0144  Y = Y.';
0145  Nsum = Nsum.';
0146 end
0147  
0148 %% % Recursive moving average code before Jos trick:
0149 % W = zeros(N,M);
0150 % Z = X(1:Wwidth,:);
0151 % inan = isnan(Z);
0152 % Z(inan) = 0;
0153 % Z = sum(Z,1);
0154 % W(F+1,:) = sum(inan,1);
0155 % Y(F+1,:) = Z;
0156 % % Recursive sum
0157 % for n = F+2:N-F
0158 %  Z = [Y(n-1,:); X(n+F,:); -X(n-F-1,:)];
0159 %  inan = isnan(Z);
0160 %  Z(isnan(Z)) = 0;
0161 %  Y(n,:) = sum(Z,1);
0162 %  W(n,:) = W(n-1,:) + inan(2,:) - inan(3,:);
0163 % end
0164 % W = Wwidth - Nnan;
0165 % Y(W==0) = NaN;
0166 % Y = Y./Wwidth;
0167 
0168 %% Remove first and last (these values are left intact above and can be
0169 %% highly erroneous)
0170 if 1
0171     Y(1) = nan;
0172     Y(end) = nan;
0173 end
0174

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005