欠損サンプルのあるデータ セットのピリオドグラム

1610 年の冬の間、ガリレオ・ガリレイは木星で最大の 4 つの衛星の動きを観察しました。ガリレオは天候が許す限り衛星の位置を記録しました。この観察記録を使用して衛星の 1 つであるカリストの軌道周期を推定します。

カリストの角度位置は分角で測定されています。曇天などの気象条件によって記録できなかったデータは、NaN を使用して指定されます。最初の観察日は 1 月 15 日です。観測時間の datetime 配列を生成します。

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

plomb を使用してデータのパワー スペクトルを推定します。オーバーサンプリング係数を 10 に指定します。結果の周波数を 1 日の秒数の逆数で表します。

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

findpeaks を使用してスペクトルのピークで目立つもののみについて位置を特定します。パワー スペクトルをプロットして、ピークを表示します。

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

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

カリストの軌道周期 (日単位) を最大エネルギーの周波数の逆数として求めます。その結果と NASA の公表値との差は 1% 未満です。

Period = 1/f0

Period = 
16.6454

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

PercentDiscrep = 
-0.2613

参考

findpeaks | plomb