function [meanout, errout] = mtmean(timeseries, nw)

% Andrew Sornborger: implements Thomsons multitaper mean for a single timeseries
%   Version: 1.0 (04/24/2006)
% mtmean: multitaper harmonic analysis mean of a timeseries using 
%  time-bandwidth 'nw'.

% length of timeseries
len = size(timeseries,1);
% number of tapers
n = 2 * nw - 1;
% replicate timeseries for tapering
totxtex = repmat(timeseries, [1 n]);
% calculate slepian tapers
seq = dpss(len(1), nw, n, 'calc');
% calculate averages of slepians
seqavg = sum(seq, 1);
jm = sum(totxtex .* seq);
% tapered mean estimate
mn = sum(jm .* seqavg)/sum(seqavg.^2);
% error of mean estimate
errout = 1/(2*nw-3) * sum((jm - mn * seqavg) .* conj(jm - mn * seqavg));
meanout = mn;
