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
    tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
    fldOut=fldOut+tmp1;
    tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
    fldCnt=fldCnt+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	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	%
11	%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	tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
56	fldOut=fldOut+tmp1;
57	tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
58	fldCnt=fldCnt+tmp1;
59	end
60
61	fldCnt(fldCnt==0)=NaN;
62	fldOut=fldOut./fldCnt;
63
64	if ~doCycle;
65	fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:);
66	%consistent with old version bug (one point offset)
67	% fldOut=fldOut(halfWindow:end-halfWindow,:,:,:);
68	end;
69
70	%switch dimensions order back to original order
71	fldOut=permute(fldOut,perm2to1);
72
73	%switch back to gcmfaces format if needed
74	if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;
75
76
77