不等間隔サンプル信号のリサンプリングおよびフィルター処理

ライブスクリプトを開く

ある人は、うるう年の 2012 年にポンド単位で体重を記録しました。この人は自分の体重を毎日は記録していないため、そのデータは不等間隔です。信号アナライザー アプリを使用して、記録された体重を前処理して調べます。アプリによって、等間隔グリッドに信号を内挿することによって欠損データ点を埋められます (この手順によって、信号に小さいギャップがある場合にのみ最適な結果が得られます)。

データを読み込み、測定値をキログラムに変換します。データファイルには、欠損している読み取りが NaN に設定されています。27 個のデータ点が欠損しています。そのほとんどは 8 月の 2 週間のストレッチの間です。

wt = datetime(2012,1,1:366)';

load weight2012.dat
wgt = weight2012(:,2)/2.20462;

validpoints = ~isnan(wgt);
missing = wt(~validpoints);
missing(15:26)

ans = 12x1 datetime
   09-Aug-2012
   10-Aug-2012
   11-Aug-2012
   12-Aug-2012
   15-Aug-2012
   16-Aug-2012
   17-Aug-2012
   18-Aug-2012
   19-Aug-2012
   20-Aug-2012
   22-Aug-2012
   23-Aug-2012

MATLAB® timetable にデータを保存します。欠損点を削除します。変動に集中できるよう、DC 値を削除します。時間情報を配列 duration に変換するために、最初の時間点を減算します。詳細については、信号アナライザーでサポートされるデータ型を参照してください。

wgt = wgt(validpoints);
wgt = wgt - mean(wgt);

wt = wt(validpoints);
wt = wt - wt(1);

wg = timetable(wt,wgt);

"信号アナライザー" を開いて timetable をディスプレイにドラッグします。[表示] タブで、[スペクトル] をクリックしてスペクトルビューを開きます。[時間] タブの [マーカーの表示] を選択します。[時間範囲] を 200 日と 250 日に設定して、欠損しているストレッチを拡大します。

信号テーブル内の信号を右クリックし、Duplicate を選択します。信号テーブルの信号を右クリックし、[名前の変更] を選択して、コピーの名前を Preprocessed に変更します。信号テーブル内で前処理済み信号を選択し、[アナライザー] タブの [前処理] をクリックします。[関数] ギャラリーで、[リサンプリング] を選択します。表示された [関数パラメーター] パネルで、以下のパラメーターを指定します。

リサンプリング法 — Sample Rate
周波数単位 — cycles/day
サンプルレート — 1
内挿法 — Shape Preserving Cubic 法

[適用] をクリックしてから [すべて確定] をクリックし、結果を保存して前処理モードを終了します。その名前の隣にあるチェックボックスを選択して表示領域上でリサンプリングされた信号を重ね合わせます。

1 年全体のデータを表示するためにズームアウトします。[スペクトル] タブで、漏れを最大値に設定します。元の信号およびリサンプリングされた信号のスペクトルはほとんどの周波数にかなり一致します。スペクトルは、2 つの目立ったピークを示します。1 つは 0.14 サイクル/日あたりで、もう 1 つは非常に低い周波数です。ピークの配置を良くするには、[表示] タブの [データカーソル] をクリックし、Two を選択します。ピーク上にカーソルを配置します。各カーソルの周波数フィールドに置き、その位置のより正確な値を取得します。

中周波数のピークは 0.143 = 1/7 サイクル/日にあり、週単位のサイクルに相当します。
低周波数のピークは 0.005 サイクル/日にあり、210 日のサイクルに相当します。

カーソルを削除するには、[データカーソル] をクリックします。表示から元の信号を除去します。Preprocessed 信号をフィルター処理して、サイクルの影響を除去します。信号テーブル内で前処理済み信号が選択された状態で、[アナライザー] タブの [前処理] をクリックします。前処理モード内で、次を行います。

低周波数のサイクルを除去するために信号をハイパスフィルター処理します。[関数] ギャラリーから [ハイパス] を選択します。表示された [関数パラメーター] パネルで、0.05 cycles/day の通過帯域周波数を入力します。他のパラメーターの既定値を使用します。[適用] をクリックします。
週次周期を除去するために信号をバンドストップフィルター処理します。[関数] ギャラリーから [バンドストップ] を選択します。[関数パラメーター] パネルで、0.135 cycles/day の低域通過帯域周波数と 0.15 cycles/day の高域通過帯域周波数を入力します。他のパラメーターの既定値を使用します。[適用] をクリックします。

結果を検査して [すべて確定] をクリックします。両方の信号を表示にプロットします。

前処理済みの信号は、元の信号より変動が少なく表示されます。信号の形状は、人の体重が冬季より夏季の月の方が少なく変動しますが、リサンプリングのアーティファクトである可能性があることを示しています。Preprocessed 信号の信号テーブルのエントリで [情報] 列のアイコンをクリックして、実行される前処理ステップを参照します。

選択した設定すべてを含む、前処理ステップの概要をすべて参照するには、[アナライザー] タブの [関数の生成] をクリックします。生成された関数が MATLAB® エディターで開きます。

function [y,ty] = preprocess(x,tx)
%  Preprocess input x
%    This function expects an input vector x and a vector of time values
%    tx. tx is a numeric vector in units of seconds.
%    Follow the timetable documentation (type 'doc timetable' in
%    command line) to learn how to index into a table variable and its time
%    values so that you can pass them into this function.

% Generated by MATLAB(R) 9.13 and Signal Processing Toolbox 9.0.
% Generated on: 03-Jan-2022 15:47:27

targetSampleRate = 1.1574074074074073e-05;
[y,ty] = resample(x,tx,targetSampleRate,'pchip');

Fs = 1/mean(diff(ty)); % Average sample rate
y = highpass(y,5.787e-07,Fs,'Steepnes',0.85,'StopbandAttenuation',60);

y = bandstop(y,[1.5625e-06 1.73611111111111e-06],Fs,'Steepness',0.85,'StopbandAttenuation',60);
end

不等間隔サンプル信号のリサンプリングおよびフィルター処理

参考

アプリ

関数

関連する例

詳細