% Calculating the best period estimate of entropy data based on Harmonic Analysis and 
% its 95% Confidence Interval based the 1000 bootstrap entropy samples
s=load('entropy_bootstrap.txt');
[m,n] = size(s); 
ysize=n;
for bootstrapIndex=1:1000
entropydata=s(bootstrapIndex,:);
y=entropydata';
periodindex=0;
for period=3.30:0.001:3.90   % This Iteration may run for a long time. The increment and the range is adjustable
periodindex=periodindex+1;
for residue=1:ysize
     xsin(residue,periodindex)= sin(2*pi/period*residue);
     xcos(residue,periodindex)= cos(2*pi/period*residue);
end
xsintem= xsin(:,periodindex);
xcostem= xcos(:,periodindex);
onestem=ones(ysize,1);
X=[onestem xsintem xcostem];
[b,bint,r,rint,stats] = regress(y,X);
Rsquare(periodindex)=stats(1);
end
[xlocation,ylocation]=max(Rsquare); 
periodestimator=3.65+(ylocation-1)*0.001;  % The period of the harmonic fitting curve with the largest R-square statistics is reported  
period_estimator(bootstrapIndex,1)=periodestimator; 
end
standarderror= std(period_estimator)/(1000)^0.5;
CIlow=mean(period_estimator)-1.96*standarderror;
CIhigh=mean(period_estimator)+1.96*standarderror;