누락된 샘플이 있는 데이터 세트의 주기도

갈릴레이 갈릴레오는 1610년 겨울 동안 목성의 가장 큰 위성 4개의 움직임을 관측했습니다. 날씨가 허락하는 경우 갈릴레오는 위성의 위치를 기록했습니다. 갈릴레오의 관측을 토대로 위성 중 하나인 칼리스토의 궤도 주기를 추정합니다.

칼리스토의 각위치는 분각(minutes of arc) 단위로 측정되었습니다. 구름 낀 날씨로 인해 누락된 데이터는 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