function [f,P,Pred_norm,Pwhite_norm,P95,P95chi2]=myspectrum(time,data)

dt = abs(mean(diff(time))); % This calculates the average sample spacing
N = length(data); 
X = fft(data)/N;
Y = X(2:ceil(length(X)/2)); %Take the first half of the fft, excluding the zero frequency.
n = length(Y);
C2 = Y.*conj(Y); % shorthand for C2 = A.^2 + B.^2 where A = real(Y) and B = imag(Y)

f = 0:(0.5/(n)):0.5'; % frequencies, from 0 to the Nyquist frequency
f = f(2:end); % get rid of the zero frequency which gives you infinite power
samf = 1/dt; % sampling frequency is one over the sample spacing
f = f * samf; % samf is the sampling frequency
P = C2/f(1); % This is the power density. The first frequency (f(1)) is the same thing as the spacing of frequencies.

r = corrcoef(data(1:end-1),data(2:end)); % Determine the red noise memory time (autolag-1 correlation coefficient)
a = r(2);
% a = (r(2)+sqrt(r(2)))/2;  % Alternative, more precise estimate of (a).

Pred = (1-a^2)./(1-2*a.*cos(2*pi.*f*dt)+a.^2);
scale = sum(P) / sum(Pred);  % This scales the red noise spectrum to the data spectrum.
Pred_norm = scale*Pred;

% Confidence levels
P95 = -log(.025)*Pred_norm;  % Note that the sqrt term in D.B.'s notes is actually for >>2 degrees of freedom,
                             % as it assumes a Gaussian distibution.  See
                             % discussion of this in the Muller and
                             % MacDonald chapter.

Pwhite_norm = mean(Pred_norm)*ones(n,1); % A flat line (on a log log plot) with the same mean power as the red noise and the data.

dof = 10; % Edit this line to adjust the degrees of freedom.
P95chi2 = chi2inv(.975,dof)/dof*Pred_norm; % This is the standard way to do the calculation for a guassian distribution.


