function theSpec = specav(theSpectra, ...
								theFrequencies, theFreqWayPoints, ...
								theTimeSpan, nPieces, theConfidenceLevel)
								
% specav -- Spectral averaging.
%  theSpec = specav(theSpectra, theFrequencies, ...
%   theFreqWayPoints, theTimeSpan, nPieces, theConfidenceLevel) averages
%   theSpectra, a Matlab "spectrum" array. whose rows correspond to
%   theFrequencies.  The bands are delimited by theFreqWayPoints.  When
%   theTimeSpan (time-wise piece-length) is non-zero, the auto and cross
%   spectra are scaled to be "variance preserving".  The effective
%   number-of-pieces defaults to 1.  The result is returned in a
%   single Matlab struct.
%
%  Note: A Matlab "spectrum" array has eight-columns, containing the
%   following: [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc].  We use just the
%   first five and recompute all the statistics.  See "help spectrum".
 
% Copyright (C) 1999 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 14-Jan-1999 07:41:29.
% Updated    25-Jan-1999 17:14:47.
% Updated    02-Feb-1999 10:15:55.

if nargin < 3, help(mfilename), return, end

if nargin < 4, theTimeSpan = 0; end
if nargin < 5, nPieces = 1; end
if nargin < 6, theConfidenceLevel = 0.95; end

if ischar(theSpectra)
	theSpectra = evalin('caller', theSpectra);
end

if ischar(theFrequencies)
	theFrequencies = evalin('caller', theFrequencies);
end

if ischar(theFreqWayPoints)
	theFreqWayPoints = evalin('caller', theFreqWayPoints);
end

if ischar(theTimeSpan)
	theTimeSpan = evalin('caller', theTimeSpan);
end

if ischar(nPieces)
	nPieces = evalin('caller', nPieces);
end

if ischar(theConfidenceLevel)
	theConfidenceLevel = evalin('caller', theConfidenceLevel);
end

% Rearrange if simple auto-spectrum.

s = theSpectra;
[m, n] = size(s);

if n == 2
	temp = zeros(m, 8);
	temp(:, [1 6]) = s(:, [1 2]);
	s = temp;
	[m, n] = size(s);
end

% Frequencies.

f = theFrequencies(:);
fi = theFreqWayPoints(:);

df = diff(fi);   % Band-widths.
nf = df ./ mean(diff(f));     % Number of frequencies per band.
cf = 0.5 * (fi(1:end-1) + fi(2:end));   % Central frequencies.

% Average the spectra.

as = diff(interp1(f, cumtrapz(f, s(:, 1:5)), fi));
as = as ./ (df * ones(1, 5));

% Extract individual spectra.

Pxx = as(:, 1); Pyy = as(:, 2); Pxy = as(:, 3);
Txy = as(:, 4); Cxy = as(:, 5);

% Recompute the statistics, following the description
%  of Marlene A. Noble, Spectral confidence levels,
%  in: Mary Hunt, PROSPECT, a sofware package for
%  spectral analysis, Woods Hole Oceanographic
%  Institution, date unknown.

cdeg = nPieces.* nf;
adeg = 2*cdeg;
degf = 2*cdeg - 2;
degf(degf < 1) = 1;
alpha = 1 - theConfidenceLevel;

% Auto-spectral confidence factors.

bmax = adeg ./ mdchi(alpha/2, adeg);
bmin = adeg ./ mdchi(1 - alpha/2, adeg);

% Coherence confidence level and phase-angle confidence.

factor = alpha;
cohlevel = sqrt(1 - factor.^(1 ./ (cdeg - 1)));
gamma2 = abs(Cxy).^2;

% Missing parenthesis in the notes regarding (1-gamma2).
%  But the example seems okay.

stat = mdsti(factor, degf);
num = 1-gamma2;
den = degf.*gamma2;
arg = stat .* sqrt(num./den);
arg(abs(arg) > 1) = NaN;
dphik = asin(arg) .* 180 ./ pi;   % Degrees.

% Transfer-function confidence interval: x-to-y and y-to-x.

dtf1 = sqrt((mdfi(theConfidenceLevel, 2, degf) .* ...
				(1-gamma2) .* Pyy) ./ (Pxx .* (cdeg-1)));

dtf2 = sqrt((mdfi(theConfidenceLevel, 2, degf) .* ...
				(1-gamma2) .* Pxx) ./ (Pyy .* (cdeg-1)));

% Variance-preserving spectra if non-zero time-span.

if theTimeSpan ~= 0
	as(:, 1:3) = as(:, 1:3) .* (cf * ones(1, 3)) .* theTimeSpan;
end

% Rename the spectra.

Pxx = as(:, 1); Pyy = as(:, 2); Pxy = as(:, 3);
Pclo = bmin; Pchi = bmax;   % Correct order?
Txyc1 = dtf1; Txyc2 = dtf2;
Cxyc = cohlevel; Cxycp = dphik;

% Pack the spectra and confidences into a struct.

ss = [];
ss.Pxx = [Pxx Pclo Pchi];   % Pxx, Pxxc-lower, Pxxc-upper.
ss.Pyy = [Pyy Pclo Pchi];   % Pyy, Pyyc-lower, Pyyc-upper.
ss.Pxy = [Pxy Pclo Pchi];   % Pxy, Pxyc-lower, Pxyc-upper.
ss.Txy = [Txy Txyc1 Txyc2];   % Txy, Txyc1, Tyxc1.
ss.Cxy = [Cxy Cxyc Cxycp];   % Coherence, coh-level, phase-uncertainty.
ss.DOF = adeg;   % Degrees-of-freedom for spectra.
ss.Freq = [cf df];   % Central frequencies and bandwidths.

% NaNs where not finite.

ss.Pxx(~isfinite(ss.Pxx)) = NaN;
ss.Pyy(~isfinite(ss.Pyy)) = NaN;
ss.Pxy(~isfinite(ss.Pxy)) = NaN;
ss.Txy(~isfinite(ss.Txy)) = NaN;
ss.Cxy(~isfinite(ss.Cxy)) = NaN;
ss.DOF(~isfinite(ss.DOF)) = NaN;
ss.Freq(~isfinite(ss.Freq)) = NaN;

if nargout > 0
	theSpec = ss;
else
	assignin('caller', 'ans', ss)
	disp(ss)
end
