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')