%% Exercise 3: Milankovitch Cycles in ODP 659 Benthic Oxygen Isotope Record
% Load the data
clear
data = load('odp659data.txt');

%%
% Plotting the data
plot(data(:,1),data(:,2))
xlabel('Alter (kyr)'),ylabel('d18O (permille vs. PDB)'),title('ODP 659')

%%
% Create three segments of the time series, 0-1,000 kyr, 1,000-2,700 kyr
% and >2,700 kyr
part1 = data(1:251,:);
part2 = data(252:669,:);
part3 = data(670:1194,:);

%%
% Plotting the segments
subplot(3,1,1),plot(part1(:,1),part1(:,2))
subplot(3,1,2),plot(part2(:,1),part2(:,2))
subplot(3,1,3),plot(part3(:,1),part3(:,2))

%%
% Min and max values of the time vector 1.960 to 5,268.000 kyrs
min(data(:,1))
max(data(:,1))

%%
% Mean sampling interval
intv1 = diff(part1(:,1));
mean(intv1)
clf, plot(intv1)
xlabel('Alter (kyr)'),ylabel('delta-Alter (kyr)'),title('ODP 659')

%%
% We choose a sampling interval of 4.4 kyrs, which is close to the mean
% sampling interval and does not increase or decrease the number of data
% points
range(data(:,1))/length(data(:,1))

%%
% We create a new time vector
t = 2 : 4.4 : 5200;
t = t';

%%
% Interpolating the data
data_L(:,1)=t;
data_L(:,2) = interp1(data(:,1),data(:,2),t,'linear');

%%
% Create segments again
part1_L=data_L(1:227,:);
part2_L=data_L(228:614,:);
part3_L=data_L(614:1182,:);

%%
% Plot the results
clf
subplot(3,1,1),plot(part1_L(:,1),part1_L(:,2))
subplot(3,1,2),plot(part2_L(:,1),part2_L(:,2))
subplot(3,1,3),plot(part3_L(:,1),part3_L(:,2))

%%
% Calculate means and variances of the time series, save in the variable
% STATISTICS
statistics(1,1)=mean(part1_L(:,2));statistics(1,2)=var(part1_L(:,2));
statistics(2,1)=mean(part2_L(:,2));statistics(2,2)=var(part2_L(:,2));
statistics(3,1)=mean(part3_L(:,2));statistics(3,2)=var(part3_L(:,2));
statistics

%%
% Detrending the data
part1_L(:,2)=detrend(part1_L(:,2));
part2_L(:,2)=detrend(part2_L(:,2));
part3_L(:,2)=detrend(part3_L(:,2));

%% 
% Calculate periodogram of segment no. 1
clf
[Pxx,f] = periodogram(part1_L(:,2),[],1024,1/4.4);
magnitude = abs(Pxx);
plot(f,magnitude);
xlabel('Frequency')
ylabel('Power')
title('Power Spectral Density Estimate 0 - 1,000 kyrs')

%%
% Calculate periodogram of segment no. 2
[Pxx,f] = periodogram(part2_L(:,2),[],1024,1/4.4);
magnitude = abs(Pxx);
plot(f,magnitude);
xlabel('Frequency')
ylabel('Power')
title('Power Spectral Density Estimate 1,000 - 2,700 kyrs')

%%
% Calculate periodogram of segment no. 3
[Pxx,f] = periodogram(part3_L(:,2),[],1024,1/4.4);
magnitude = abs(Pxx);
plot(f,magnitude);
xlabel('Frequency')
ylabel('Power')
title('Power Spectral Density Estimate 2,700 - 5,200 kyrs')

%% 
% Detrending the full linearly-interpolated time series
data_L(:,2) = detrend(data_L(:,2));

%%
% Evolutionary autostrectrum
spectrogram(data_L(:,2),128,90,1024,1/4.4);

%%
% Full Autospectrum
[Pxx1,f1] = periodogram(data_L(:,2),[],1024,1/4.4);
plot(f1,abs(Pxx1))

%%
% Designing bandstop filters to remove the 100, 40 and 20 kyr cycles
[b1,a1] = butter(4,[0.005 0.015]/0.1136,'stop');
[h1,w1] = freqz(b1,a1,512);
f1= (1/4.4)*w1/(2*pi);
magnitude1 = abs(h1);

[b2,a2] = butter(5,[0.015 0.035]/0.1136,'stop');
[h2,w2] = freqz(b2,a2,512);
f2= (1/4.4)*w2/(2*pi);
magnitude2 = abs(h2);

[b3,a3] = butter(7,[0.04 0.06]/0.1136,'stop');
[h3,w3] = freqz(b3,a3,1024);
f3= (1/4.4)*w3/(2*pi);
magnitude3 = abs(h3);

%%
% Filter characteristics
subplot(1,3,1),plot(f1,magnitude1),grid,axis([0 0.1 0 1])
xlabel('Frequency'), ylabel('Magnitude')
title('Bandstop Filter 100 kyr')

subplot(1,3,2),plot(f2,magnitude2),grid,axis([0 0.1 0 1])
xlabel('Frequency'), ylabel('Magnitude')
title('Bandstop Filter 40 kyr')

subplot(1,3,3),plot(f3,magnitude3),grid,axis([0 0.1 0 1])
xlabel('Frequency'), ylabel('Magnitude')
title('Bandstop Filter 20 kyr')

%%
% Zero-phase forward and reverse filtering
data_F(:,1) = t;
data_F(:,2) = filtfilt(b1,a1,data_L(:,2));
data_F(:,3) = filtfilt(b2,a2,data_L(:,2));
data_F(:,4) = filtfilt(b3,a3,data_L(:,2));

%%
% Filter output in the time domain
subplot(3,1,1),plot(data_F(:,1),data_F(:,2))
xlabel('Alter (kyr)'),ylabel('d18O (permille vs. PDB)'),
title('ODP 659 ohne 100 kyr Zyklus')

subplot(3,1,2),plot(data_F(:,1),data_F(:,3))
xlabel('Alter (kyr)'),ylabel('d18O (permille vs. PDB)'),
title('ODP 659 ohne 40 kyr Zyklus')

subplot(3,1,3),plot(data_F(:,1),data_F(:,4))
xlabel('Alter (kyr)'),ylabel('d18O (permille vs. PDB)'),
title('ODP 659 ohne 20 kyr Zyklus')

%%
% Filter output in the frequency domain
[Pxx1,f1] = periodogram(data_F(:,2),[],1024,1/4.4);
subplot(1,3,1),plot(f1,abs(Pxx1))
xlabel('Frequency'), ylabel('Power')
title('ODP 659 ohne 100 kyr Zyklus')

[Pxx2,f2] = periodogram(data_F(:,3),[],1024,1/4.4);
subplot(1,3,2),plot(f2,abs(Pxx2))
xlabel('Frequency'), ylabel('Power')
title('ODP 659 ohne 40 kyr Zyklus')

[Pxx3,f3] = periodogram(data_F(:,4),[],1024,1/4.4);
subplot(1,3,3),plot(f3,abs(Pxx3))
xlabel('Frequency'), ylabel('Power')
title('ODP 659 ohne 20 kyr Zyklus')