www.gusucode.com > signal 案例源码程序 matlab代码 > signal/PeriodogramOfADataSetWithMissingSamplesExample.m

    %% Periodogram of Data Set with Irregular Sampling
% Galileo Galilei discovered Jupiter's four largest satellites in January
% of 1610 and recorded their locations every clear night until March of
% that year. Use Galileo's data to find the orbital period of Callisto, the
% outermost of the four satellites.

% Copyright 2015 The MathWorks, Inc.


%%
% Galileo's observations of Callisto's angular position are in minutes of
% arc. There are several gaps due to cloudy conditions.

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]';

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]';

plot(t,yg,'o')
xlabel('Day')

%%
% Use |plomb| to compute the periodogram of the data. Estimate the power
% spectrum up to a frequency of $0.5\;{\rm day}^{-1}$. Specify an
% oversampling factor of 10. Choose the standard Lomb-Scargle
% normalization.

[pxx,f] = plomb(yg,t,0.5,10,'normalized');

%%
% The periodogram has one clear maximum. Name the peak frequency $f_0$.
% Plot the periodogram and annotate the peak.

[pmax,lmax] = max(pxx);
f0 = f(lmax);

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

%%
% Use linear least squares to fit to the data a function of the form
%
% $$y(t)=A+B\cos2\pi f_0t+C\sin2\pi f_0t.$$
%
% The fitting parameters are the amplitudes _A_, _B_, and _C_.

ft = 2*pi*f0*t;

ABC = [ones(size(ft)) cos(ft) sin(ft)] \ yg

%%
% Use the fitting parameters to construct the fitting function on a
% 200-point interval. Plot the data and overlay the fit.

x = linspace(t(1),t(end),200)';
fx = 2*pi*f0*x;

y = [ones(size(fx)) cos(fx) sin(fx)] * ABC;

plot(t,yg,'o',x,y)
xlabel('Day')