clear all; close all;
% This replicates Figure 2 of MBH99 preprint
% Jean S/5.10.2006

residualsAD1000=load('MBH99_AD1000_res.txt');
residualsAD1820=load('MBH99_AD1820_res.txt');

% Notice that the "time-bandwidth product" is 2
% (i.e., the number of spheroidal sequences is 
% then 5. One can experiment with this parameter
% to see how "robust" Mann's conclusions are if
% the uncertainty variances indeed are modified
% according to the spectral estimates.

PxxAD1000=pmtm(residualsAD1000,2);
PxxAD1820=pmtm(residualsAD1820,2);

% The first and the last coefficient seem
% not to be in the spectrum figure. This has something
% to do with the implementation of the MT method,
% I suppose. So I just ignore them below.

years=1902:1980;
l=length(PxxAD1000)-2;
frequency=[0:l-1]/(2*(l-1));
% "white noise variance" (i.e., mean power level)
S0_AD1000=mean(PxxAD1000(2:end-1));
S0_AD1820=mean(PxxAD1820(2:end-1));

figure; 
subplot(2,1,1);
plot(years,residualsAD1000,'r'); hold on;
plot(years,residualsAD1820,'b');
axis([1900 1980 -0.8 0.6]);
legend('AD1000','AD1820',0);
grid on;

subplot(2,1,2);
semilogy(frequency,PxxAD1000(2:end-1)/S0_AD1000,'r'); hold on;
semilogy(frequency,PxxAD1820(2:end-1)/S0_AD1820,'b');
axis([0 0.5 0 6]);
xlabel('FREQUENCY (cycle/yr)');
ylabel('RELATIVE VARIANCE');
legend('AD1000','AD1820',0);
grid on;

print -dpng figure2.png
save 'MBH99_AD1000_spectrum.txt' PxxAD1000 -ASCII 
save 'MBH99_AD1820_spectrum.txt' PxxAD1820 -ASCII
