function y = nanmovmean( x, mlag, mspan )

% (C) Demetrio Ferro demetrio.ferro@upf.edu
% Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International Public License

% the function applies a sliding nanmean to the first dimension of input x
% centered at each mlag and with width msize, allowing overlap for msize > mlag

% e.g. size(x,1)=12,                      mlag=2       msize=6
% x : [ 1 2 3 4 5 6 7 8 9 10 11 12 ]     [ - - ]   [ - - - - - - ] 
% 
% y(1)=nanmean([ 1  2  3  4                        ]);        
% y(2)=nanmean([ 1  2  3  4  5  6                  ]);
% y(3)=nanmean([       3  4  5  6  7  8            ]);
% y(4)=nanmean([             5  6  7  8  9 10      ]);
% y(5)=nanmean([                   7  8  9 10 11 12]);
% y(6)=nanmean([                         9 10 11 12]);

% choose msize and mlag as multiple of 2 or fix this code accordingly

if mod(mlag,2)==1 
    error('the parameter mlag has to be divisible by 2.')
end
if isscalar(mspan)
if mod(mspan,2)==1
    error('the parameter mspan has to be divisible by 2.')
else
    mspan_bf=mspan/2;
    mspan_af=mspan/2;
end
elseif isvector(mspan) && length(mspan)==2
    mspan_bf=mspan(1);
    mspan_af=mspan(2);
else
    error('specify a valid mspan value');
end

if isvector(x)
    L=length(x);
    y=nan(1,ceil(L/mlag));
    for ll=1:ceil(L/mlag)
        ctsmp=max(((ll-1)*mlag)-mspan_bf+1+mlag/2,1) : min(((ll)*mlag)+mspan_af-mlag/2,L);
        y(ll)=nanmean(x(ctsmp));
    end
elseif ismatrix(x)
    [L,M]=size(x);
    y=nan(ceil(L/mlag),M);
    for ll=1:ceil(L/mlag)
        ctsmp=max(((ll-1)*mlag)-mspan_bf+1+mlag/2,1) : min(((ll)*mlag)+mspan_af-mlag/2,L); 
        y(ll,:)=nanmean(x(ctsmp,:),1);
    end
end
end