大阪大学医学部 Python会

Now is better than never.

ZapLineを用いたLFPに対する電源由来のアーティファクト除去 (Python)

2022-05-21(Sat) - Posted by 山本 in 技術ブログ    tag:Python tag:Neuroscience

Contents

    ZapLineとは

    ZapLineはEEGやMEGなどのデータにおける電力系統 (power line)に由来するアーティファクトを除去するアルゴリズムである．アーティファクトの周波数を指定し（残念ながら指定が必要），アーティファクト成分のみを除去する時空間的フィルタを求めることを目的とする．やや手法が込み入っていたので詳細は論文のmethodを読んでほしい．

    ZapLineの論文：

    Cheveigné, Alain de. 2020. “ZapLine: A Simple and Effective Method to Remove Power Line Artifacts.” NeuroImage.

    meegkitでZapLineを使用する

    ZapLineはMATLABでの実装が著者により公開されている．

    この記事ではmeegkitというPythonライブラリにおける実装を使用する．documentとgithubへのリンクは次の通り：

    必要なライブラリのinstall

    pipにより，pyedflib, pyriemann, meegkitをinstallする．

    pip install pyedflib
pip install pyriemann
pip install git+https://github.com/nbara/python-meegkit.git

    pyedflibは.bdfファイルを読み込むために使用する．pyriemannはmeegkitで使用されており，別途installが必要となった．

    ライブラリのimport

    以下ではmeegkitのexampleを参考にしてZapLineの使用法を見ていく．まずライブラリをimportする．

    In [1]:
    import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
plt.rcParams.update({'axes.spines.top': False, 'axes.spines.right': False})

from meegkit import dss
from meegkit.utils import unfold
from meegkit.utils.denoise import demean

import pyedflib

    データのdownload

    解析に使用するデータをdownloadする．著者実装で使用されていたEEG_Cat_Study4_Resting_S1.bdfをRaw Empirical EEG Data - Trujillo et al. (2017) Frontiers in NeuroscienceのページでDownloadする．

    あるいはexampleのようにdss_line_data.npyを使用してもよい．

    In [2]:
    file = pyedflib.EdfReader("EEG_Cat_Study4_Resting_S1.bdf")

    bdfファイルを読み込む．メソッドの詳細はpyedflibを参照．

    In [3]:
    n = file.signals_in_file # 73
signal_labels = file.getSignalLabels() # len(signal_labels): 73

sigbufs = np.zeros((n, file.getNSamples()[0])) # sigbufs.shape: (73, 130560)
for i in np.arange(n):
    sigbufs[i, :] = file.readSignal(i)

#header = file.getSignalHeaders()
sfreq = file.getSampleFrequencies()[0] # sample frequency

    データの一部だけをsliceし，demean関数でcenteringする．

    In [4]:
    X = sigbufs[:67, :20000]
X = demean(X.T)

    パワースペクトル密度 (power spectral density; PSD)をscipy.signal.welchで推定する．

    In [5]:
    freq, Pxx = signal.welch(unfold(X), sfreq, nperseg=500, axis=0, 
                         return_onesided=True)
    In [6]:
    fig, ax = plt.subplots(1, 2, figsize=(8, 4),
                       sharey=True, constrained_layout=True)
ax[0].semilogy(freq, Pxx, lw=.5)
ax[0].set_title('All channels')
ax[0].set_xlabel('Frequency [Hz]')
ax[0].set_ylabel(r'PSD [$V^2$/Hz]')
ax[0].grid(True, which="major")
ax[0].grid(True, which="minor", linestyle='dashed')
ax[0].set_xlim(freq[0], freq[-1])

ax[1].semilogy(freq, Pxx.mean(1), lw=2, label="original")
ax[1].set_title('Mean')
ax[1].set_xlabel('Frequency [Hz]')
ax[1].grid(True, which="major")
ax[1].grid(True, which="minor", linestyle='dashed')
ax[1].set_xlim(freq[0], freq[-1])
ax[1].legend()
plt.show()

    左が使用したチャネルのPSDであり，右はその平均である．60Hz付近に電源に由来するアーティファクトが見える．これがアーティファクトと言える理由は60Hz丁度にあり，ピーク幅が狭いなどが挙げられる．これを除く方法を次節で紹介する．

    ZapLineの実行

    交流電力の周波数は50Hz or 60Hzなのでいずれかを指定してZapLineを実行する．もちろんこれ以外の周波数も除去できる．

    In [7]:
    fline = 60 # line frequency

# Apply dss_line (ZapLine)
out, _ = dss.dss_line(X, fline, sfreq, nremove=3, nfft=1024)

    Power of components removed by DSS: 0.00

    nremoveは除去アルゴリズムを実行する回数であり，収束するまで実行するdss.dss_line_iter関数もある．ただ収束しない場合もあるので，そうした場合にはdss.dss_lineを使用する．

    結果を描画する．

    In [8]:
    freq, Pxx_out = signal.welch(unfold(out), sfreq, nperseg=500, axis=0, 
                             return_onesided=True)
    In [9]:
    fig, ax = plt.subplots(1, 2, figsize=(8, 4), 
                       sharey=True, constrained_layout=True)
ax[0].semilogy(freq, Pxx_out, lw=.5)
ax[0].set_title('All channels')
ax[0].set_xlabel('Frequency [Hz]')
ax[0].set_ylabel(r'PSD [$V^2$/Hz]')
ax[0].grid(True, which="major")
ax[0].grid(True, which="minor", linestyle='dashed')
ax[0].set_xlim(freq[0], freq[-1])

ax[1].semilogy(freq, Pxx_out.mean(1), lw=2, label="clean")
ax[1].semilogy(freq, np.abs(Pxx - Pxx_out).mean(1), lw=1, 
               label="removed")
ax[1].set_title('Mean')
ax[1].set_xlabel('Frequency [Hz]')
ax[1].grid(True, which="major")
ax[1].grid(True, which="minor", linestyle='dashed')
ax[1].set_xlim(freq[0], freq[-1])
ax[1].legend()
plt.show()

    左がZapLineを適応した全チャネルのPSDであり，右は青がその平均，オレンジは元データとの差分の絶対値である．一つ上の図と見比べると，概ねアーティファクト成分のみを除けていることがわかる．

    Tags Python  Neuroscience