Periodogram of Data Set with Missing Samples

Open Live Script

Galileo Galilei observed the motion of Jupiter's four largest satellites during the winter of 1610. When the weather allowed, Galileo recorded the satellites' locations. Use his observations to estimate the orbital period of one of the satellites, Callisto.

Callisto's angular position is measured in minutes of arc. Missing data due to cloudy conditions are specified using NaNs. The first observation is dated January 15. Generate a datetime array of observation times.

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

obsv = datetime(1610,1,14+(1:length(yg)));

plot(yg,'o')

ax = gca;
nights = [1 18 32 46];
ax.XTick = nights;
ax.XTickLabel = char(obsv(nights));
grid

Figure contains an axes object. The axes contains a line object which displays its values using only markers.

Estimate the power spectrum of the data using plomb. Specify an oversampling factor of 10. Express the resulting frequencies in inverse days.

[pxx,f] = plomb(yg,obsv,[],10,'power');
f = f*86400;

Use findpeaks to determine the location of the only prominent peak of the spectrum. Plot the power spectrum and show the peak.

[pk,f0] = findpeaks(pxx,f,'MinPeakHeight',10);

plot(f,pxx,f0,pk,'o')
xlabel('Frequency (day^{-1})')
title('Power Spectrum and Prominent Peak')
grid

Figure contains an axes object. The axes object with title Power Spectrum and Prominent Peak, xlabel Frequency (day toThePowerOf - 1 baseline ) contains 2 objects of type line. One or more of the lines displays its values using only markers

Determine Callisto's orbital period (in days) as the inverse of the frequency of maximum energy. The result differs by less than 1% from the value published by NASA.

Period = 1/f0
Period = 
16.6454

NASA = 16.6890184;
PercentDiscrep = (Period-NASA)/NASA*100
PercentDiscrep = 
-0.2613

See Also

|