Обработка сигнала с отсутствующими выборками
Рассмотрим вес человека как зарегистрированный (в фунтах) в течение високосного 2012 года. Человек не записывал свой вес каждый день. Вы хотели бы изучить периодичность сигнала, но прежде чем вы сможете сделать это, вы должны позаботиться о недостающих данных.
Загрузите данные и преобразуйте измерения в килограммы. Пропущенные показания устанавливаются в NaN. Определите, сколько точек отсутствует.
load('weight2012.dat') wgt = weight2012(:,2)/2.20462; daynum = 1:length(wgt); missing = isnan(wgt); fprintf('Missing %d samples of %d\n',sum(missing),max(daynum))
Missing 27 samples of 366
Присвоение значений отсутствующим точкам с помощью resample. По умолчанию resample делает оценки с помощью линейной интерполяции. Постройте график исходных и интерполированных показаний. Масштабирование в дни с 200 по 250, которые содержат около половины недостающих точек.
wgt_orig = wgt; wgt = resample(wgt,daynum); plot(daynum,wgt_orig,'.-',daynum,wgt,'o') xlabel('Day') ylabel('Weight (kg)') axis([200 250 73 77]) legend('Original','Interpolated') grid
Определите, является ли сигнал периодическим, проанализировав его в частотный диапазон. Найдите длительность цикла, измерив время в неделях. Вычесть среднее значение, чтобы сконцентрироваться на колебаниях.
Fs = 7;
[p,f] = pwelch(wgt-mean(wgt),[],[],[],Fs);
plot(f,p)
xlabel('Frequency (week^{-1})')
grid
Заметьте, как вес человека колеблется еженедельно. Есть ли заметный шаблон от недели до недели? Исключить последние два дня года, чтобы получить 52 недели. Переупорядочить измерения в соответствии с днем недели.
wgd = reshape(wgt(1:7*52),[7 52]); plot(wgd') xlabel('Week') ylabel('Weight (kg)') dweek = datetime([repmat([2012 1],7,1) (1:7)'],'Format','eeee'); legend(string(dweek),'Location','northwest') grid
Сглаживайте колебания с помощью фильтра, который подходит для полиномов низкого порядка к подмножествам данных. В частности, установите его в соответствии с кубическими полиномами к наборам семи дней.
wgs = sgolayfilt(wgd',3,7); plot(wgs) xlabel('Week') ylabel('Weight (kg)') title('Smoothed Weight Fluctuations') legend(string(dweek),'Location','southeast'); grid
Этот человек, как правило, ест больше, и, таким образом, весит больше, в выходные. Проверьте путем вычисления ежедневных средств.
for jk = 1:7 fprintf('%s mean: %5.1f kg\n',dweek(jk),mean(wgd(jk,:))) end
Sunday mean: 76.2 kg Monday mean: 75.7 kg Tuesday mean: 75.2 kg Wednesday mean: 74.9 kg Thursday mean: 75.1 kg Friday mean: 75.3 kg Saturday mean: 75.8 kg