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
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