matlab_class/gcmfaces_misc/runmean.m

function [fldOut]=runmean(fldIn,halfWindow,dim,varargin);
%object:    compute running mean window ('rmw') over a dimension
%input:     fldIn is the field to which the rmw will be applied
%           halfWindow is the half width of the rmw
%           dim is the dimension over which the rmw will be applied
%optional:  doCycle states whether the boundary condition is cyclic (1) or
%               not (0; default). If doCycle==0, the no. of averaged points
%               decreases from 1+2*halfWindow to halfWindow at the edges.
%output:    fldOut is the resulting field
%
%notes:     - NaNs are discarded in the rmw, implying that an average is
%           computed if the rmw contains at least 1 valid point.
%           - setting halfWindow to 0 implies fldOut=fldIn.

%determine and check a couple things
fld_isa_gcmfaces=isa(fldIn,'gcmfaces');
if fld_isa_gcmfaces;
    if dim<3;
        error('for gcmfaces objects runmean excludes dim=1 or 2');
    end;
end;

if nargin==3; doCycle=0; else; doCycle=varargin{1}; end;

%switch to array format if needed
if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;

%switch dim to first dimension
sizeIn=size(fldIn); 

perm1to2=[1:length(sizeIn)]; 
perm1to2=[dim perm1to2(find(perm1to2~=dim))];
perm2to1=[[1:dim-1]+1 1 [dim+1:length(sizeIn)]]; 
sizeIn2=sizeIn(perm1to2);

fldIn=permute(fldIn,perm1to2); 

sizeCur=size(fldIn);
if ~doCycle;
    %add NaNs at both edges
    sizeCur(1)=sizeCur(1)+2*halfWindow;
    fldIn2=NaN*ones(sizeCur);
    fldIn2(halfWindow+1:end-halfWindow,:,:,:)=fldIn;
    fldIn=fldIn2; clear fldIn2;
end;

%create mask and remove NaNs:
fldMsk=~isnan(fldIn);
fldIn(isnan(fldIn))=0;
fldCnt=0*fldIn;

%apply the running mean
fldOut=zeros(sizeCur);
for tcur=-halfWindow:halfWindow
    %To have halfWindow*2 coeffs rather than halfWindow*2+1, centered to the current 
    %point, it is convenient to reduce the weight of the farthest points to 1/2.
    %This became necessary to get proper annual means, from monthly data, with halfWindow=6.
    if abs(tcur)==halfWindow; fac=1/2; else; fac=1; end;
    tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
    fldOut=fldOut+fac*tmp1;
    tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
    fldCnt=fldCnt+fac*tmp1;
end

fldCnt(fldCnt==0)=NaN;
fldOut=fldOut./fldCnt;

if ~doCycle;
    fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:);
    %consistent with old version bug (one point offset)
    %     fldOut=fldOut(halfWindow:end-halfWindow,:,:,:);
end;

%switch dimensions order back to original order
fldOut=permute(fldOut,perm2to1);

%switch back to gcmfaces format if needed
if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;


1	gforget	1.2	function [fldOut]=runmean(fldIn,halfWindow,dim,varargin);
2			%object: compute running mean window ('rmw') over a dimension
3			%input: fldIn is the field to which the rmw will be applied
4			% halfWindow is the half width of the rmw
5			% dim is the dimension over which the rmw will be applied
6			%optional: doCycle states whether the boundary condition is cyclic (1) or
7			% not (0; default). If doCycle==0, the no. of averaged points
8			% decreases from 1+2*halfWindow to halfWindow at the edges.
9			%output: fldOut is the resulting field
10	gforget	1.1	%
11	gforget	1.2	%notes: - NaNs are discarded in the rmw, implying that an average is
12			% computed if the rmw contains at least 1 valid point.
13			% - setting halfWindow to 0 implies fldOut=fldIn.
14
15			%determine and check a couple things
16			fld_isa_gcmfaces=isa(fldIn,'gcmfaces');
17			if fld_isa_gcmfaces;
18			if dim<3;
19			error('for gcmfaces objects runmean excludes dim=1 or 2');
20			end;
21			end;
22
23			if nargin==3; doCycle=0; else; doCycle=varargin{1}; end;
24
25			%switch to array format if needed
26			if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;
27
28			%switch dim to first dimension
29			sizeIn=size(fldIn);
30
31			perm1to2=[1:length(sizeIn)];
32			perm1to2=[dim perm1to2(find(perm1to2~=dim))];
33			perm2to1=[[1:dim-1]+1 1 [dim+1:length(sizeIn)]];
34			sizeIn2=sizeIn(perm1to2);
35
36			fldIn=permute(fldIn,perm1to2);
37
38			sizeCur=size(fldIn);
39			if ~doCycle;
40			%add NaNs at both edges
41			sizeCur(1)=sizeCur(1)+2*halfWindow;
42			fldIn2=NaN*ones(sizeCur);
43			fldIn2(halfWindow+1:end-halfWindow,:,:,:)=fldIn;
44			fldIn=fldIn2; clear fldIn2;
45			end;
46
47			%create mask and remove NaNs:
48			fldMsk=~isnan(fldIn);
49			fldIn(isnan(fldIn))=0;
50			fldCnt=0*fldIn;
51
52			%apply the running mean
53			fldOut=zeros(sizeCur);
54			for tcur=-halfWindow:halfWindow
55	gforget	1.3	%To have halfWindow2 coeffs rather than halfWindow2+1, centered to the current
56			%point, it is convenient to reduce the weight of the farthest points to 1/2.
57			%This became necessary to get proper annual means, from monthly data, with halfWindow=6.
58			if abs(tcur)==halfWindow; fac=1/2; else; fac=1; end;
59	gforget	1.2	tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
60	gforget	1.3	fldOut=fldOut+fac*tmp1;
61	gforget	1.2	tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
62	gforget	1.3	fldCnt=fldCnt+fac*tmp1;
63	gforget	1.2	end
64	gforget	1.1
65	gforget	1.2	fldCnt(fldCnt==0)=NaN;
66			fldOut=fldOut./fldCnt;
67	gforget	1.1
68	gforget	1.2	if ~doCycle;
69			fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:);
70			%consistent with old version bug (one point offset)
71			% fldOut=fldOut(halfWindow:end-halfWindow,:,:,:);
72			end;
73	gforget	1.1
74	gforget	1.2	%switch dimensions order back to original order
75			fldOut=permute(fldOut,perm2to1);
76	gforget	1.1
77	gforget	1.2	%switch back to gcmfaces format if needed
78			if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;
79	gforget	1.1
80
81