8. 周波数解析

8.1 パワースペクトラムを見る

R_median (61 files)をFiles to processにドラッグし、RUNを押してFrequency→Power spectrum density (Welch)を選択します。

Edit...を押します。

Matlab’s FFT defaultsとSave average PSD valuesを選択し、OKを押します。

Runを押します。

PSD: 1/4000ms Avg,Power (MEG)を右クリックしてPower spectrumを選択します。

パワースペクトラムが表示されます。magnetometerとplanar gradiometerで値が違うので２種類表示されています。西日本なので60Hzの整数倍の雑音が観察されています。対数表示を改めるにはPower、Magnitudeなどを選択します。

とりあえずこのデータは要らないのでPSD: 1/4000ms Avg,Power (MEG)を右クリックし、File→Deleteを選択します。

はい(Y)を押します。

8.2 加算平均波形の時間周波数解析

正式な時間周波数解析は各epoch毎に計算した後、平均したものをいいますが、とりあえず簡易版です。

Avg: R_median (61 files)をFile to processにドラッグし、RUNを押します。

Frequency→Hilbert transformを選択します。

Runを押します。

Magnitudeを右クリックして2D Layout (maps)を選択します。

タブをDisplayにしてPowerを選択したのち、2D Layout画像を右クリックしてColormap: Timefreq→Maximum: Globalを選択します。

どこかのチャンネルをクリックすると、クリックしたチャンネルの時間周波数マップが表示されます。

チャンネルのウィンドウを閉じ、2D Layout画像を右クリックしTopographyを選択します。

Norm Gradを押します。この画面は最初に定義すると、次からは出てこないようですが、出る時もあります・・・。

トポグラフィが表示されます。時間を22ms、周波数帯域をbetaにします。

トポグラフィ画面を右クリックしてMaximum: Localを選択します。

Magnitude (MEG)を右クリックし、3D Sensor capを選択します。

3D表示されます。

周波数成分を細かく見ていきます。

ウィンドウをすべて閉じ、Avg: R_median (61 files)をFiles to processにドラッグし、RUNを押してFrequency→Time-frequency (Morlet wavelets)を選択します。

Edit...を押します。

OKを押します。

Runを押します。

Power, 1:1:60Hz (MEG)を右クリックし2D Layout (maps)を選択します。

2D Layoutを右クリックしてMaximum: Globalを選択します。

Hide edge effectsにチェックを入れます。

時間の両端、低周波数帯域が除去された時間周波数マップになります。

22ms、20HzにしてTopographyを選択します。

Norm Gradを押します。

Topographyが表示されました。1Hz刻みで細かく周波数帯域を観察できるようになりました。

8.3 各エポックごとの時間周波数マップの加算平均作成

ウィンドウは全部閉じます。

エポックごと切り出されたR median (61 files)をFiles to processにドラッグし、Time-frequency (Morlet wavelets)を選択します。

Edit...を押します。

Save averaged...を選択し、OKを押します。

Runを押します。計算処理に時間がかかります。

Hide edge effectsのチェックを外し、Avg,Power, 1:1:60Hz (MEG)を右クリックして2D Layout (maps)を選択します。

画面を右クリックしMaximum: Custom...を選択します。

Maximumを0.5にしてOKを押します。

Topographyを選択します。

時間を22.0ms、周波数を30Hzとしました。

ウィンドウをすべて閉じ、Avg,Power,1:1:60Hz(MEG)をFiles to processにドラッグし、RUNを押してStandardize→Z-score(static)を選択します。

Runを押します。

Avg,Power,1:1:60Hz(MEG)|zscoreを2D Layout (maps)表示します。

チャンネルをクリックし22ms、20Hzにしました。赤が事象関連同期、青が事象関連脱同期です。但し、誘発磁場成分は除去していません。

チャンネルの赤十字線をクリックして39.9ms、49Hzの時点で3D Sensor capを作成しました。

8.4 各エポックから電流源推定した結果の時間周波数マップを加算平均

R_median (61 files)をFile to processにドラッグしてSources Compute sourcesを選択し、Runを押します。

Sourceを選択します。

Frequency→Time-frequency (Morlet wavelets)を選択します。

Use scouts time seriesのチェックを外し、Edit...を押します。

Group in frequency bands (Hz)を選択し、OKを押し、Runを押します。計算時価はそこそこかかります。

Avg,Power,Freqbandsを右クリックしDisplay on cortexを選択します。

22ms、alpha 8-12Hzを選択してMaximum Custom..を選択し、Maximum:Customeの最大値を8にします。

任意の位置で右クリックし、Source: Time-frequency Shiftを選択します。

指定した点の時間周波数マップが表示されます。

Avg,Power,FreqBandsをFile to processにドラッグしてRUNを押し、Standardize→Event related pertubation (ERS/ERD)を選択し、Runを押します。

Avg,Power,FreqBands|ersdを右クリックしDisplay on cortexを選択します。

22ms、alpha帯域、Colormapのminimumとmaximumを-10～10としてみました。

gamma1帯域にして任意の点を選択し、その部位の時間周波数マップを表示させました。

画面を右クリックしColormap:Timefreq→Colormapでrbwを選択しました。

カラースケールの色が変わりました。