load Huybers07data.txt
figure
time = Huybers07data(:,1)
Benthic = Huybers07data(:,4);
% Split the originald data into five bins, but use only the first half (up to 1,000,000 years)

k = find(time<1000);
time = time(k);
Benthic = Benthic(k);
n = length(Benthic);
binsize = floor(n/5);

bin1 = Benthic(1:binsize);
bin2 = Benthic(binsize+1:binsize*2);
bin3 = Benthic(binsize*2+1:binsize*3);
bin4 = Benthic(binsize*3+1:binsize*4);
bin5 = Benthic(binsize*4+1:binsize*5);

time = time(1:length(bin1));

[f1,P1,Pred_norm1,Pwhite_norm1,P951,P95chi21]=myspectrum(time,bin1);
[f2,P2,Pred_norm2,Pwhite_norm2,P952,P95chi22]=myspectrum(time,bin2);
[f3,P3,Pred_norm3,Pwhite_norm3,P953,P95chi23]=myspectrum(time,bin3);
[f4,P4,Pred_norm4,Pwhite_norm4,P954,P95chi24]=myspectrum(time,bin4);
[f5,P5,Pred_norm5,Pwhite_norm5,P955,P95chi25]=myspectrum(time,bin5);

f = mean([f1;f2;f3;f4;f5]);
P = mean([P1';P2';P3';P4';P5']);
Pred_norm = mean([Pred_norm1;Pred_norm2;Pred_norm3;Pred_norm4;Pred_norm5]);
P95chi = mean([P95chi21;P95chi22;P95chi23;P95chi24;P95chi25]);
subplot(1,3,1);
loglog(f,P);
hold on
loglog(f,Pred_norm,'r');
loglog(f,P95chi,'g');

% do this again but use only the second part (past 1,000,000 years)

clear;
load Huybers07data.txt

time = Huybers07data(:,1)
Benthic = Huybers07data(:,4);
k = find(time>1000);
time = time(k);
Benthic = Benthic(k);
n = length(Benthic);
binsize = floor(n/5);

bin1 = Benthic(1:binsize);
bin2 = Benthic(binsize+1:binsize*2);
bin3 = Benthic(binsize*2+1:binsize*3);
bin4 = Benthic(binsize*3+1:binsize*4);
bin5 = Benthic(binsize*4+1:binsize*5);

time = time(1:length(bin1));

[f1,P1,Pred_norm1,Pwhite_norm1,P951,P95chi21]=myspectrum(time,bin1);
[f2,P2,Pred_norm2,Pwhite_norm2,P952,P95chi22]=myspectrum(time,bin2);
[f3,P3,Pred_norm3,Pwhite_norm3,P953,P95chi23]=myspectrum(time,bin3);
[f4,P4,Pred_norm4,Pwhite_norm4,P954,P95chi24]=myspectrum(time,bin4);
[f5,P5,Pred_norm5,Pwhite_norm5,P955,P95chi25]=myspectrum(time,bin5);

f = mean([f1;f2;f3;f4;f5]);
P = mean([P1';P2';P3';P4';P5']);
Pred_norm = mean([Pred_norm1;Pred_norm2;Pred_norm3;Pred_norm4;Pred_norm5]);
P95chi = mean([P95chi21;P95chi22;P95chi23;P95chi24;P95chi25]);

subplot(1,3,2);
loglog(f,P);
hold on
loglog(f,Pred_norm,'r');
loglog(f,P95chi,'g');

% Now bin the spectral data.



clear;
load Huybers07data.txt

time = Huybers07data(:,1)
Benthic = Huybers07data(:,4);

[f,P,Pred_norm,Pwhite_norm,P95,P95chi2]=myspectrum(time,Benthic);

for i = 1:floor(length(f)/5)
    j = (i-1)*5+1;
    k = j+4;
    fm(i)=mean(f(j:k)');
    Pm(i)=mean(P(j:k));
    Pred_normm(i)=mean(Pred_norm(j:k));
    P95chi2m(i)=mean(P95chi2(j:k));
end

subplot(1,3,3);
loglog(fm,Pm);
hold on
loglog(fm,Pred_normm,'r');
loglog(fm,P95chi2m,'g');


