function [f, A, B] = myfft_coef(data, SPD)

% sine and cosine coefficients from a time series

% OUTPUTS

% f = the frequency vector- discrete frequencies resolved
% A = cosine coefficient for the component with frequency f and assuming the time series starts at time zero  
% B= sine coefficient, as above
%        
% the time domain would be given by

t=0:1/SPD: (length(data)-1)*(1/SPD);


% 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);
A=real(coef); %cosine
A(NF)=A(NF)/2; %afjust the nyquist frequency
B=-imag(coef); %sine, not sure why the negative is needed



