% This program is used to compute the fft power spectrum for a real signal
% format s=fftspectrum(signal)
% Note: signal should be given as a column vector.
%  
% .............................................................................
% References: 
%
% Zhou YH, SalsteinDA, Chen JL, Revised atmospheric excitation function series 
%     related to Earth's variable rotation under consideration of surface topography, 
%     Journal of Geophysical Research, 111, D12108, doi:10.1029/2005JD006608, 2006.
%
% Jenkins GM, Watts DG, Spectrum analysis and its applications. Holden-Day
%     Press, California, 1968.
% .............................................................................
%
% Yonghong Zhou, Feb.2009, at Shanghai Astronomical Observatory 
%
function s=fftspectrum(signal)
%--- Power Spectral Analysis (using FFT) --
%We use FFT to get the estimation of Power Spectral Density.
%The units on the y-axis are u.^2/(unit of frequency), if 
%the signal datax has units of u associated with it. 
%
n0=length(signal);
n1=2^(nextpow2(n0));
Yn=fft(signal,n1);
%
Pyy=Yn.*conj(Yn)/n0;
%
%Because we do (n1 points) FFT, n should equal n1;
n=length(Yn);
%
interval=input('sampling interval = ');
nyquist=1/(2*interval);
%
nn=n/2;
ff=[0:nn]/nn*nyquist;
ff=ff';
PP(1)=Pyy(1);
PP(2:nn)=2*Pyy(2:nn);
PP(nn+1)=Pyy(nn+1);
PP=PP';
PP=10*log10(PP);
s=[ff PP];
figure,plot(s(:,1),s(:,2));
xlabel('frequency');
% For example, if the unit of sampling interval is "day",
% then, the unit of frequency is "cycles per day". 
%
ylabel('db'); 
% fft power spectral density estimates in decibels