不等間隔サンプル信号のリサンプリングおよびフィルター処理
ある人は、うるう年の 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
参考
アプリ
関数
関連する例
- 相関する信号間の遅延の検出
- ウィンドウの漏れを変化させることでトーンを分解する
- パーシステンス スペクトルを使用した干渉の検出
- 複素包絡線を使用した変調と復調
- 再代入したスペクトログラムを使用したリッジの検出と追跡
- 音楽信号からの音声の抽出
- 独自の関数を使用した飽和信号のクリップ除去
- 振動信号の包絡線スペクトルの計算
- クジラの歌からの関心領域の抽出