function [freq, pow]=f_powerspec(SDATA, SF)

%* use Matlab's fft function to do a simple spectral analysis of isotope
% time-series; DFT
% DC component of x is fftx(1) 
% fftx(1+nfft/2)> is Nyquist frequency component
%
% SDATA: data value at corresponding sample times
% SF:    sample frequency
% freq:  corresponding frequency vector
% 

spad= 2^(nextpow2(length(SDATA)));   % figure closest power of 2 so fft can pad
spec=fft(SDATA,spad);          % fourier transform
N=ceil((spad+1)/2);            % number of unique points 
spec=spec(1:N);                % result of fft is redundant, take first half plus one to get DC
pow=abs(spec)';                % magnitudes, returns the complex modulus where needed
pow=pow/length(SDATA);         % scale the fft so no dependence on length(SDATA)
pow=pow.^2;
pow=pow*2;                     % fix for half thrown away
pow(1)=pow(1)/2;               % no fix for DC component, it is unique
if ~rem(spad,2)                % no fix for Nyquist if it is there
    pow(end) = pow(end)/2; 
end
freq=(0:N-1)*SF/spad;             % frequency vector