Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

欠損サンプルのある信号の処理

うるう年の 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

該当者の体重が週次でどのように変動したかに注目します。週単位で注目すべきパターンがあるでしょうか。年の最後の 2 日間は 52 週にするため消去します。曜日別に測定値を並べ替えます。

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

データのサブセットを低次多項式に近似するフィルターを使用して、変動を滑らかにします。特に、7 日間のセットを 3 次多項式に近似するように設定します。

wgs = sgolayfilt(wgd',3,7);

plot(wgs)
xlabel('Week')
ylabel('Weight (kg)')
title('Smoothed Weight Fluctuations')

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

該当者は、週末に他の曜日よりも多く食べるため体重も増加する傾向があります。日次平均値を計算して確認します。

for jk = 1:7
    fprintf('%3s mean: %5.1f kg\n', ...
        datestr(datenum(2012,1,jk),'ddd')',mean(wgd(jk,:)))
end
Sun mean:  76.2 kg
Mon mean:  75.7 kg
Tue mean:  75.2 kg
Wed mean:  74.9 kg
Thu mean:  75.1 kg
Fri mean:  75.3 kg
Sat mean:  75.8 kg

参考

|

関連するトピック