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