function [S,f] = powspec(x,varargin)

% POWSPEC
%   For a discrete time-series X, returns the Power Spectrum
%   of the sequence.
% AUTHOR
%   (c) S.Angus, Jan 2007 (UNSW)
% INPUTS
%   (X) the time-series (vector) to analyse.
% OPTIONS
%   (...,S) optionally, the resultant spectra can be smoothed
%     using the moving average of S pieces of data.
%   (...,X1X2) an OLS can be performed and graphed if needed
%     on a selection of the data, represented by a two-value
%     vector X1X2 for the starting and finishing data to
%     consider.
% OUTPUTS
%   [S,f] The spectra S and the corresponding frequencies F.
% See also FFT SMOOTH

% --------------------------------------------------------------- %
s = 0;
x1x2 = 0;
% --------------------------------------------------------------- %

if nargin > 1, s = varargin{1}; end
if nargin > 2, x1x2 = varargin{2}; end

n = length(x);
F = fft(x);				% discrete Fourier Transform
cut = floor((n-1)/2)+1;	% only take LHS of series
G = F(2:cut);			% omit constant-shift term f(0)
S = [G.*conj(G)/n^2]';	% get and normalise magnitude of complex compenents
f = [2:cut]'/n;			% get frequencies
Ss = S;

if s, Ss = smooth(S,s); end

% plot
figure(gcf),clf
loglog(f,Ss,'k-'),grid on, axis square,hold on
xlabel 'f'
ylabel 'S(f)'
title(sprintf('Power Spectrum (smoothed,%d pts)',s))
pretty_plot(gcf,gca)
disp(sprintf('Points: %d',cut-1));

if x1x2
	% prepare data for OLS
	x = log(f(x1x2(1):x1x2(2)));
	y = log(Ss(x1x2(1):x1x2(2)));
	res = ols(y,[ones(length(x),1) x]);
	loglog(f(x1x2(1):x1x2(2)),exp(res.yhat),'r-')
	disp(sprintf('Slope: %6.4f',res.beta(2)));
end


