% This program is used to compute the multitaper power spectrum for a real signal
% format s=mtspectrum(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.
%
% Percival D. B. and A. T. Walden, Spectral Analysis for Physical Applications: 
%     Multitaper and Conventional Univariate Techniques. Cambridge, England: Cambridge 
%     University Press£¬1993.
% .............................................................................
%
% Yonghong Zhou, Feb 2009, at Shanghai Astronomical Observatory 
%
function s=mtspectrum(signal)
s1=mtspec(signal);
n=length(s1);
interval=input('sampling interval = ');
nyquist=1/(2*interval);
cc=[0:(n-1)]/n*nyquist;
cc=cc';
s=[cc s1];
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'); 
% Multitaper power spectral density estimates in decibels






% This program computes the 4pi-multitapered power spectrum 
% of (real-valued) time series.
%
function s1=mtspec(signal)
n=length(signal);
m_taper=7;
m=m_taper;
disp('fft spectrum...')
if m==0
   s=fft(signal);
   s=s.*conj(s)/n;
   n=ceil(n/2);
   s1=10*log10(s(1:n));
else
   s=zeros(n,1);
   ss=s;
   nw=(m+1)/2;  % m=2*nw-1
   [v,lambda]=dpss(n,nw);
   for j=1:m
       ss=fft(v(:,j).*signal);
       ss=ss.*conj(ss)/lambda(j);
       s=s+ss;
   end
   n=ceil(n/2);
   s1=10*log10(s(1:n)/m);
end