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

    %% Detect Periodicity in a Signal with Missing Samples
% Consider the weight of a person as recorded (in pounds) during the leap
% year 2012. The person did not record their weight every day. You would
% like to study the periodicity of the signal, even though some data points
% are missing.

% Copyright 2015 The MathWorks, Inc.


%%
% Load the data and convert the measurements to kilograms. Missed readings
% are set to |NaN|. Determine how many points are missing.

load(fullfile(matlabroot,'examples','signal','weight2012.dat'))

wgt = weight2012(:,2)/2.20462;

fprintf('Missing %d samples of %d\n',sum(isnan(wgt)),length(wgt))

%%
% Determine if the signal is periodic by analyzing it in the frequency
% domain. The Lomb-Scargle algorithm is designed to handle data with
% missing samples or data that have been sampled irregularly.
% 
% Find the cycle durations, measuring time in weeks.

[p,f] = plomb(wgt,7,'normalized');

plot(f,p)
xlabel('Frequency (week^{-1})')

%%
% Notice how the person's weight oscillates weekly. Is there a noticeable
% pattern from week to week? Eliminate the last two days of the year to get
% 52 weeks. Reorder the measurements according to the day of the week.

wgd = reshape(wgt(1:7*52),[7 52])';

plot(wgd)
xlabel('Week')
ylabel('Weight (kg)')

q = legend(datestr(datenum(2012,1,1:7),'dddd'));
q.Location = 'NorthWest';

%%
% Smooth out the fluctuations using a filter that fits low-order
% polynomials to subsets of the data. Specifically, set it to fit cubic
% polynomials to sets of seven days.

wgs = sgolayfilt(wgd,3,7);

plot(wgs)
xlabel('Week')
ylabel('Smoothed weight (kg)')

q = legend(datestr(datenum(2012,1,1:7),'dddd'));
q.Location = 'SouthEast';

%%
% This person tends to eat more, and thus weigh more, during the weekend.
% Verify by computing the daily means. Exclude the missing values from the
% calculation.

for jk = 1:7
    wgm = find(~isnan(wgd(:,jk)));
    fprintf('%3s mean: %5.1f kg\n', ...
        datestr(datenum(2012,1,jk),'ddd')',mean(wgd(wgm,jk)))
end