function [fs, pows] = myfft_s(data, SPD,smoothnumber)

% display the power spectrum of discrete modes after being exposed to a
% one tenth cosine bell and smoothed in the spectral domain

% 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) 
%smoothnumber = number of spectral estimates to smooth over


%make the input data a column vector
[NN, N]=size(data);
if NN==1
    data=data'; %transpose of the vector is a row vector
end


N=length(data); %number of points in domain
t=1:N;
t=t-1; %make the time domain start at zero
T=length(t);
otcb=ones(T,1);


% figure(19); clf;
% plot(t, data, 'g', t, cosb'.*data, 'r')

%do the fft
Yfft=fft(otcb.*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
variancereductionofwindow=var(otcb.*data)/var(data)

for n=1:rounddown((NF-1)/smoothnumber);
    fs(n)=mean(f(2+(n-1)*smoothnumber:1+n*smoothnumber));
    pows(n)=sum(pow(2+(n-1)*smoothnumber:1+n*smoothnumber));
end

if rounddown((NF-1)/smoothnumber)*smoothnumber+1 < NF
    fnum=length(rounddown((NF-1)/smoothnumber)*smoothnumber+2:NF); %number of spectral estimates in the final band
    %warning(['the last spectral bin has only ', num2str(fnum),  '  spectral estimates in it'])
    fs(n+1)= mean(f(rounddown((NF-1)/smoothnumber)*smoothnumber+2:NF));
    pows(n+1)=sum(pow(rounddown((NF-1)/smoothnumber)*smoothnumber+2:NF));
end

 %figure(1); clf;
 %plot(fs,pows, 'b*', f, pows, 'b')
 %xlabel('FREQUENCY (CYCLES PER YEAR)', 'FontSize', 18)
 %ylabel('SPECTRAL POWER (sum of discrete modes is the variance)', 'FontSize', 18)


