function [I,kF,kS]=tperiod(x,fig)
%        COMPUTING PERIODOGRAM AND DETECTING SIGNIFICANT PERIODICITIES
%                      USING FISHER'S AND SIEGEL'S TEST
%        Additive model: x(t)=y(t)+e(t) is assumed where
%        y(t) ... deterministic component with periodicities to be detected
%        e(t) ... error component, e: IID(0,sigma^2)
% =============================================================================
% function [I,kF,kS]=tperiod(x,fig)
% INPUT:
% x   ... is an n x 1 or 1 x n data vector of observations
% fig ... figure handle for the peridogram with emphasized significant
%         harmonics
%         optional: []=default=no plot
% OUTPUT:
% I ... periodogram: n2 x 1 or 1 x n2 vector sized according to x
%                    n2=fix((n-1)/2)
% kF ... indices of significant harmonics using Fisher's test
% kS ... indices of significant harmonics using Siegel's test

% (C) 1997-1999, V. Vesely, Masaryk University of Brno (Czech Republic)
% Rev. 1, Dec. 22, 1999
%
error(nargchk(1,2,nargin));
if nargin<2, fig= []; end

if min(size(x))>1, error('x is not a vector'); end
n=length(x);
m=fix((n-1)./2) ;             % length of the periodogram
if m<5, error('x must be longer than 10'); end
m=min(m,1000);
I=fft(x);
I=(2./n)*abs(I(2:m+1)).^2;    % Computing the periodogram
II=I;
lF=[]; lS=[]; k=[];
while m>4                     % loop to detect significant periodicities
   [gf,ts]=fisher(m);         % Fisher's and Siegel's 95% percentile
   Y=II./sum(II);
   [W,kmax]=max(Y);           % Fisher's statistic
   T=sum(max(0,Y-0.6*gf));    % Siegel's statistic
   if isempty(lF) & W<=gf, lF=length(k); end
   if isempty(lS) & T<=ts, lS=length(k); end
   if W>gf | T>ts
      k=[k,kmax(1)];
      II(kmax(1))=0; m=m-1;   % scratch component just detected
   else
      break;
   end
end
if ~isempty(fig)
   stem(I,'--');
   keephold=ishold;
   hold on; stem(k,I(k)); stem(k,I(k),'*');
   if ~keephold, hold off; end
   if nargout>2 & lS>lF
      title('PERIODOGRAM (Siegel''s test)');
   else
      title('PERIODOGRAM (Fisher''s test)');
   end
end
kF = k(1:lF);
kS = k(1:lS);
