function [f, pow, phase] = myfft(data, SPD)

% display the power spectrum of discrete modes

% OUTPUTS

% f = the frequency vector- discrete frequencies resolved
% pow = spectral power defined such that the sum of the discrete modes is the variance in the original data 
% phase = phase of the signal at each frequency
%               to reconstruct the data, add the phase to the cosine of the 
%               root of the power PLUS the phase  (a sin in the data = cos(t-pi/2)) 


% INPUTS

% data = a vector of the data on an evenly spaced domain
% SPD = samples per day, or other unit of time you want the frequencies in
%            terms of (reciprocal of) 


N=length(data); %number of points in domain

%do the fft
Yfft=fft(data);
f=SPD*(0:1:N/2)/N; %definition of frequencies, should be multiplied by samples per day
NF=length(f);
coef(1:NF)=Yfft(1:NF)/(N/2);
pow(1:NF)=((abs(Yfft(1:NF)/(N/2))).^2)/2; %POWER PER UNIT FREQUENCY
pow(NF)=((abs(Yfft(NF)/(N/2))).^2)/4; %the nyquist frequency must be adjusted
phase(1:NF)=angle(coef); %defined such that the original data is cos(t + phase) --- i.e. the sin would have phase -pi/2
figure(1); clf;
plot(f(2:end),pow(2:end), 'b*', f(2:end), pow(2:end), 'b')
xlabel('FREQUENCY (CYCLES PER DAY)', 'FontSize', 18)
ylabel('SPECTRAL POWER (sum of discrete modes is the variance)', 'FontSize', 18)


