%% Resample a Nonuniformly Sampled Data Set
% Use the data recorded by Galileo Galilei in 1610 to determine the orbital period of Callisto, the outermost of Jupiter's four largest satellites.
% Galileo observed the satellites' motion for six weeks, starting on 15 January. The observations have several gaps because Jupiter was not visible on cloudy nights. Generate a datetime array of observation times.

t = [0 2 3 7 8 9 10 11 12 17 18 19 20 24 25 26 27 28 29 31 32 33 35 37 ...
    41 42 43 44 45]'+1;

yg = [10.5 11.5 10.5 -5.5 -10.0 -12.0 -11.5 -12.0 -7.5 8.5 12.5 12.5 ...
    10.5 -6.0 -11.5 -12.5 -12.5 -10.5 -6.5 2.0 8.5 10.5 13.5 10.5 -8.5 ...
    -10.5 -10.5 -10.0 -8.0]';

obsv = datetime(1610,1,15+t);

%
% Resample the data onto a regular grid using a sample rate of one observation per day. Use a moderate upsampling factor of 3 to avoid overfitting.

fs = 1;

[y,ty] = resample(yg,t,fs,3,1);

%
% Plot the data and the resampled signal.

plot(t,yg,'o',ty,y,'.-')
xlabel('Day')

%
% Repeat the procedure using spline interpolation and displaying the observation dates. Express the sample rate in inverse days.

fs = 1/86400;

[ys,tys] = resample(yg,obsv,fs,3,1,'spline');

plot(t,yg,'o')
hold on
plot(ys,'.-')
hold off

ax = gca;
ax.XTick = t(1:9:end);
ax.XTickLabel = char(obsv(1:9:end));

%
% Compute the periodogram power spectrum estimate of the uniformly spaced, linearly interpolated data. Choose a DFT length of 1024. The signal peaks at the inverse of the orbital period.

[pxx,f] = periodogram(ys,[],1024,1,'power');
[pk,i0] = max(pxx);

f0 = f(i0);
T0 = 1/f0

plot(f,pxx,f0,pk,'o')
xlabel('Frequency (day^{-1})')

% Copyright 2012 The MathWorks, Inc.