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