一様にサンプリングされていない信号のリサンプリングおよびフィルター処理

ライブスクリプトを開く

ある人は、うるう年の 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® タイムテーブルにデータを保存します。欠損点を削除します。変動に集中できるよう、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 に変更します。Preprocessed 信号を選択したままにします。[アナライザー] タブの [前処理] ▼ をクリックして、[リサンプリング] を選択します。表示される [リサンプリング] タブで、1 cycles/day のサンプルレートを入力して、Shape Preserving Cubic 方法を選択します。[リサンプリング] をクリックします。その名前の隣にあるチェックボックスを選択して表示領域上でリサンプリングされた信号を重ね合わせます。

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

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

カーソルを削除するには、[データカーソル] アイコンをクリックします。表示から元の信号を除去します。Preprocessed 信号をフィルター処理して、サイクルの影響を除去します。

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

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

選択した設定すべてを含む、前処理ステップの概要をすべて参照するには、[アナライザー] タブの [関数の生成] をクリックします。生成された関数が 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.5 and Signal Processing Toolbox 8.1.
% Generated on: 08-Jun-2018 14:35:38

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,'Steepness',0.85,'StopbandAttenuation',60);
y = bandstop(y,[1.5625e-06 1.73611111111111e-06],Fs,'Steepness',0.85,'StopbandAttenuation',60);
end