function [fstatret, sigclean, sigavg, erravg, ...
          signale, erre, signalo, erro] = ...
         harmanal(timeseries, nw, harm, pad, nstim, stt, stt2, stt3)

% function [fstatret, sigclean, sigavg, erravg,
%           signale, erre, signalo, erro] = ...
%          harmanal(timeseries, nw, harm, pad, nstim, ...
%                       'stat', 'round', 'all')
%
% Author: Andrew Sornborger (04/24/2006, Version 1.1)
%         bug fix: empty f generates a return (9/5/00)
%         enhancement: finds all significant harmonics
%                      starting with base `harm'
%         enhancement: multitaper mean subtracted and put in average estimate
%
% Implementation of multitaper harmonic analysis using F-statistic for
%   finding and removing sinusoidal discrete lines in a signal.
% This code will detect and estimate harmonics at multiples of a base frequency.
%   It will also separately extract the non-specific (average) and specific
%   (differential) dynamics (see J. Neurosci. Meth. ref. below).
%  References:
%    Biophys. J., Vol. 76 691-708 (some typos in equations, user beware)
%    Proc. IEEE, Vol. 70 1055-1096 (better, but more technical ref.)
%    NeuroImage, Vol. 18 610-621
%    J. Neurosci. Meth., Vol. 141 223-229.
%  Input:
%    timeseries: timeseries to be analyzed
%    nw: time-bandwidth product
%    harm: wavenumber of first harmonic
%    pad: number of multiples of series size to take fft of
%    nstim: number of stimuli in experiment
%    stt: 'stat' if statistical significance to be checked, 'nostat' otherwise
%    stt2: 'round' harmonics to be rounded, 'noround' otherwise
%    stt3: 'all' all harmonics to be extracted, not just periodic harmonics
%  Returns:
%    fstatret: measure of statistical significance, I usually compare this to a
%              significance level of log10(length_of_timeseries).
%    sigclean: signal cleaned of significant harmonics
%    sigavg:  estimated mean of signal
%    erravg:   estimated error of mean
%    signale:  combined significant stimulus non-specific harmonics
%    erre:     estimated error of non-specific signal
%    signalo:  combined significant stimulus specific harmonics
%    erro:     estimated error of specific signal
%  Note: sigclean + sigavg + signale + signalo = timeseries

% length of time series
len = size(timeseries, 1);

% time
t = [0:len-1];

% number of tapers
n = 2 * nw - 1;

% subtract mean from timeseries
[mtmn, mtmnerr] = mtmean(timeseries, nw);
timeseries = timeseries - mtmn;

% make copies of time series
totxtex = repmat(timeseries, [1 n]);

% generate discrete prolate spheroidal sequences
seq = dpss(len(1), nw, n, 'calc');

% calculate U_k(0)
seqavg = sum(seq, 1);

% tapered fourier transform
xf = fft(totxtex .* seq * size(seq, 1) * pad, size(seq, 1) * pad);

% calculate estimate of sinusoid amplitude
k = [1:2:n];
mu = (xf(:,k) * seqavg(k)')/sum(seqavg .* seqavg);

% calculate f-statistic for estimated sinusoid
num = (mu .* conj(mu)) * (n - 1) * sum(seqavg .* seqavg);
denom = sum((xf - mu * seqavg) .* conj(xf - mu * seqavg), 2);
fstat = num./denom;

% calculate log of probability distribution function of f-statistic
fstatret = (fcdf(fstat, 2, 2 * n - 2));
fstatret = log10(fpdf(fstat, 2, 2 * n - 2));

j = [2:size(fstatret,1)-1];
% find locations of significant peak maxima
if (size(stt,2) == 4 & stt == 'stat')
  f = find((fstatret(j) < -log10(len)) & (fstatret(j) < fstatret(j-1)) ...
                               & (fstatret(j) < fstatret(j+1)));
else
  f = find((fstatret(j) < fstatret(j-1)) ...
                               & (fstatret(j) < fstatret(j+1)));
end;

if (fstatret(1) < -log10(len))
  f = [1; f];
end;

% return if no maxima
if (isempty(f))
  sigclean = timeseries;
  signale = zeros(size([1:len]))';
  signalo = zeros(size([1:len]))';
  sigavg = zeros(size([1:len]))';
  erravg = 0;
  erre = 0;
  erro = 0;
  return;
end;

% determine all significant even harmonics
%fns4 = f(find(mod(f - 1, S * harm * pad) <= pad));
%fns5 = f(find(mod(f - 1, S * harm * pad) > S * harm * pad - pad));
fns4 = f(find(mod(f - 1, nstim * harm * pad) <= pad));
fns5 = f(find(mod(f - 1, nstim * harm * pad) > nstim * harm * pad - pad));
fns6 = [fns4; fns5];
fns7 = fns6(find(fns6 > harm * pad/2));
fns8 = fns7(find(fns7 < size(fstatret,1)/2));
fns9 = size(fstatret,1) - fns8([1:size(fns8,1)]) + 1;
fns10 = [fns8([1:size(fns8,1)]); fns9];
if ((size(stt2, 2) == 5) & (stt2 == 'round'))
  fns10 = round(fns10/(pad * harm)) * pad * harm;
end;

if (~isempty(fns6))
  fns11 = fns6(find(fns6 == 1));
else
  fns11 = [];
end;

%fs4 = f(find(mod(f - 1, harm * pad) <= pad ...
%             & ~(mod(f - 1, S * harm * pad) <= pad)));
%fs5 = f(find(mod(f - 1, harm * pad) > harm * pad - pad ...
%             & ~(mod(f - 1, S * harm * pad) > S * harm * pad - pad)));
fs4 = f(find(mod(f - 1, harm * pad) <= pad ...
             & ~(mod(f - 1, nstim * harm * pad) <= pad)));
fs5 = f(find(mod(f - 1, harm * pad) > harm * pad - pad ...
             & ~(mod(f - 1, nstim * harm * pad) > nstim * harm * pad - pad)));
fs6 = [fs4; fs5];
fs7 = fs6(find(fs6 > harm * pad/2));
fs8 = fs7(find(fs7 < size(fstatret,1)/2));
fs9 = size(fstatret,1) - fs8([1:size(fs8,1)]) + 1;
fs10 = [fs8([1:size(fs8,1)]); fs9];
if (size(stt2, 2) == 5 & stt2 == 'round')
  fs10 = round(fs10/(pad * harm)) * pad * harm;
end;

% construct estimated sinusoids on grid
este = zeros(size(fstat));
if (size(stt3,2) == 3 & stt3 == 'all')
  este(f) = mu(f);
else;
  este(fns10) = mu(fns10);
end;
ieste = real(ifft(este));
stdeve = sqrt(1/((len * pad)^2) * sum(1/nw * denom(fns10)/sum(seqavg .* seqavg)));

esto = zeros(size(fstat));
esto(fs10) = mu(fs10);
iesto = real(ifft(esto));
stdevo = sqrt(1/((len * pad)^2) * sum(1/nw * denom(fs10)/sum(seqavg .* seqavg)));

estavg = zeros(size(fstat));
estavg(fns11) = mu(fns11);
iestavg = real(ifft(estavg));
stdeva = sqrt(1/((len * pad)^2) * sum(1/nw * denom(fns11)/sum(seqavg .* seqavg)));

% separated signals and errors for return
if (size(stt3,2) == 3 & stt3 == 'all')
  sigclean = timeseries - ieste(1:len);
else
  sigclean = timeseries - ieste(1:len) - iesto(1:len) - iestavg(1:len);
end;
sigavg = mtmn * ones(size([1:len]))';
erravg = mtmnerr * ones(size([1:len]))';
signale = ieste(1:len);
if (~isempty(stdeve))
  erre = stdeve;
else
  erre = 0;
end;
signalo = iesto(1:len);
if (~isempty(stdevo))
  erro    = stdevo;
else
  erro = 0;
end;

fstatret = abs(fstatret);
