function [avg,Lh] = wrunavg(varargin)
% [avg,Lh] = wrunavg(data, L, ends, weight)
% ----------
% weighted running average of length L on data. Lh returns the inverse of
% the sum of the square of the weights (used for df calculations)
% ends tag tells how to handle the ends:
%  - 0: nothing is done, ends are handled as per default in conv2.m
%  - 1: 'periodic': the data is assume to be periodic and padded as such
%  - 2: 'antiperiodic': the data is assume to be antiperiodic and padded as such
%  - 3: as in 0 except the first and last L/2 points of avg are made equal to NaN
% weighted is a variable specifying the weighting
%  - 0: equal weights across L (i.e., unweighted average)
%  - 1: triangle weights (default)

[cax,args,nargs] = axescheck(varargin{:});
dat = args{1};
L = args{2};
ends = args{3};
if nargs == 4
    weight = args{4};
else
    weight = 1;
end

% make sure L is odd
if mod(L,2) == 0
    L = L + 1;
end;
m = (L-1)/2;

% get weighting function
l=0:m;
hp = 1 - l./m;
l=-m:-1;
hm = 1 + l./m;
h = [hm hp];

% non-weighted average if desired
if ( weight==0 ) h = ones(size(h)); end

% renormalize so that weighting yields same area under curve
h = h./sum(h);

% calculate Lh
Lh = sum(h.*h); Lh = 1./Lh;

% calculate average

if ( ends == 0 )
    avg = conv(dat, h); avg = avg((m+1):(length(avg)-m));
elseif ( ends == 3 ) % periodic
    avg = conv(dat, h); avg = avg((m+1):(length(avg)-m));
    avg(1:m) = NaN; avg(end-m+1:end) = NaN;
elseif ( ends == 1 ) % periodic
    dat1(1:m) = flipud(dat(1:m));
    dat1(m+1:m+length(dat)) = dat;
    dat1(m+length(dat)+1:m+length(dat)+m) = flipud(dat(length(dat)-m+1:length(dat)));
    avg = conv(dat1, h); avg = avg(m+m+1:length(avg)-m-m);
elseif ( ends == 2 ) % antiperiodic
    dat1(1:m) = -flipud(dat(1:m));
    dat1(m+1:m+length(dat)) = dat;
    dat1(m+length(dat)+1:m+length(dat)+m) = -flipud(dat(length(dat)-m+1:length(dat)));
    avg = conv(dat1, h); avg = avg(m+m+1:length(avg)-m-m);
end
