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

% display the power spectrum on a log scale
% the spectral power is scaled such that the area under the curve is the
% variance

% OUTPUTS

% f = the frequency vector- discrete frequencies resolved
% pow = spectral power defined such that the integral over frequency is the
%          variance of the original data set 
% 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)=(N/SPD)*((abs(Yfft(1:NF)/(N/2))).^2)/2; %POWER PER UNIT FREQUENCY
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(fignum); clf;
plot(log(f), f.*pow, 'b*', log(f), f.*pow, 'b')
num=8;
lfrequencies=log(f(2)):(log(f(NF))-log(f(2)))/(num-1):log(f(NF));
frequencies=exp(lfrequencies);
%maxperiod=1/f(2);
%minperiod=1/f(NF);
%periods=maxperiod:-(maxperiod-minperiod)/19:minperiod
periods=1./frequencies;
set(gca,'XTick',lfrequencies)
set(gca,'XTickLabel',num2str(round(100*periods')/100))
%XTick=log(1./periods);
xlabel(' PERIODS (DAYS) [on a ln(frequency axis)', 'FontSize', 18)
ylabel('SPECTRAL POWER TIMES FREQUENCY (area under curve is the power)', 'FontSize', 18)

%check to make sure that the variance in the original data equals the
%spectral variance in both the linear and log spectral space

do=2:NF; %the zeroeth mode has no variance, disclude it
VARIANCE_CHECK=[var(data);trapz(f(do),pow(do)); trapz(log(f(do)), f(do).*pow(do))]