ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで生体信号解析２】Welch法で脳波パワーを計算

脳波 EEG R言語 パワースペクトル 生体工学

前回の記事では、EDF ファイルの読み込について説明しました。今回はその続編として、「読み込んだ脳波信号を使って、Welch 法でパワースペクトルを計算し、さらに脳波バンドのパワーを求める」ということをやってみます。
　Welch 法は、脳波解析だけでなく、生体信号解析全般でよく使われる パワースペクトル（Power Spectral Density：PSD）の推定法です。コード例とともに、実際の解析でのポイントをまとめます。
　今回は 単一チャンネルの脳波解析 を想定していますので、EDF から取り出した

x    <- ch1$signal     # EEG 波形（1 チャンネル分）
fs   <- ch1$sRate      # サンプリング周波数 [Hz]
t0   <- ch1$startTime  # 計測開始時刻（POSIXct / 文字列）

をデータとして使います。

1. Welch 法とは

■ パワースペクトルを真の姿に近づける

　パワースペクトル密度（PSD）の推定では、観測された 1 本の時系列の背後に「真のパワースペクトル」が存在すると考えます。これは、その信号が従う確率過程の性質を反映したものであり、本来は各周波数成分の「平均的なパワー」を描いた、なだらかな曲線になります。
　この「平均的なパワー」という点が重要で、パワースペクトルとは理論的に、フーリエ変換した成分のパワーを「時間長で割ったうえで平均したもの」として定義されます。つまり、真のパワースペクトルは、確率過程の期待値として定義される、「平均値の曲線」です。
　ところが、実際の計算では、有限長のデータを 1 回だけ観測してFFT（高速フーリエ変換）を行い、そのパワーを図示することになります。このとき得られるペリオドグラムは、周波数に沿って大きく変動する「ギザギザした曲線」になります。
　ここで強調したいのは、このギザギザの主因が「データが有限長だから」であったり「FFT を使ったから」であったりするわけではないということです。そうではなく、各周波数でのパワーは本来「確率変数」であり、1 回の観測値は必ず上下に揺らぐという事実が、本質的な理由です。有限長であることは、その揺らぎの大きさを左右する要因にすぎません。
　PSD が本質的に「平均として定義される量」なのですから、ギザギザをならして真の姿に近づけるためには、どこかで平均をとる操作が必要になります。
　主な方法としては次の 2 つがあります。

  • 周波数方向でペリオドグラムを平滑化する方法
  • 複数のペリオドグラムを平均する方法（アンサンブル平均） の 2 種類があります。

　本来の定義にもっとも近いのは後者、すなわちアンサンブル平均です。これは「同じ条件で何度も信号を計測し、それぞれのパワースペクトルを計算し、最後に平均する」という操作です。こうすれば、確率変数として揺らぐペリオドグラムのばらつきを抑え、真のパワースペクトルの滑らかな姿を復元できます。
　しかし、生体信号では同じ条件で繰り返し計測することが難しいのが普通です。被験者の状態も環境も、時間とともに微妙に変化してしまうためです。
　そこで Welch 法では、次のような工夫が使われます。同じ確率過程に従う信号であれば、長い時系列を時間方向に分割して得られる部分列も、もとの信号と同じ統計的性質を持つはずだという考え方です。つまり、1 回しか測定していない長い信号を複数の短いセグメントに分割し、それぞれからペリオドグラムを求め、それらを平均するのです。こうして時間方向に作り出した「疑似アンサンブル」によって、真のパワースペクトルの滑らかな姿により近い推定を得ようとする方法が Welch 法なのです。

■ Welch 法の手順

　Welch 法では、長い時系列を複数の短い区間に切り分け、それぞれの区間からパワースペクトルを求め、その平均をとることで安定した推定値を得ます。具体的な流れは次のとおりです。

  1. 時系列を一定長のセグメントに分割する 例として「5 秒ごと」に区間を切り出すなど、短い部分時系列をいくつも作ります。

  2. 各セグメントに窓関数をかける 急激な端の影響を抑えるため、Hamming 窓などの窓関数を適用します。

  3. 各セグメントで FFT を行い、ペリオドグラム（パワースペクトル）を計算する

  4. すべてのセグメントのパワースペクトルを平均する 平均をとることで、確率的に揺らぐペリオドグラムのばらつきが抑えられ、滑らかなスペクトル推定が得られます。

  5. 平均したパワースペクトルを最終的な Welch 推定値として利用する

このように Welch 法は、時間方向に分割した複数のセグメントから「疑似アンサンブル」を作り、平均をとることで、真のパワースペクトルにより近い推定を得る方法です。

2. 脳波バンドの定義

　脳波の周波数バンドは、「どの周波数帯をどんな機能と結びつけて解釈するか」ということで、経験的に確立されてきた分類です。そのため、分野（睡眠・認知神経科学・臨床・動物実験など）や研究グループによって、周波数の値や名称が少しずつ異なることに注意が必要です。
　ここでは、最近のヒト脳波研究で広く使われている代表的な定義を示しつつ、バンドごとの機能的な解釈も整理しておきます。

■ 2.1 代表的なバンド定義（ヒト EEG の例）

　ヒトの頭皮 EEG では、よく次のような周波数帯が用いられます。

  • デルタ（delta）帯：おおよそ 0.5 ～ 4 Hz
  • シータ（theta）帯：おおよそ 4 ～ 8 Hz
  • アルファ（alpha）帯：おおよそ 8 ～ 12（あるいは 13）Hz
  • ベータ（beta）帯：おおよそ 13 ～ 30 Hz
  • ガンマ（gamma）帯：おおよそ 30 Hz 以上（45 Hz 前後まで、あるいは 80 Hz 以上を含めて議論する場合もある）

ここで重要なのは、「周波数帯の範囲は絶対ではない」という点です。たとえば、アルファ帯を 8–12 Hz とするか 8–13 Hz とするか、シータ帯を 4–7 Hz とするか 4–8 Hz とするかは、論文や研究者によって少しずつ異なります。
　また、高周波側については、

  • 臨床 EEG や睡眠解析では「30 Hz 以上」をひとまとめにガンマとして扱う
  • 認知神経科学の分野では「低ガンマ（low gamma: 30–60 Hz）」「高ガンマ（high gamma: 60–100 Hz）」など、さらに細かく分ける

といった違いがあります。

■ 2.2 各バンドの典型的な解釈

　以下は、ヒト EEG 研究でよく用いられる、機能との対応づけの一例です。あくまで「傾向」であり、必ずしも一対一に対応するものではありません。

  • デルタ（0.5–4 Hz 前後） 深いノンレム睡眠（徐波睡眠）で顕著に増加する、大きくてゆっくりとした波です。 集団レベルでみると、デルタパワーは睡眠の深さや、脳の「ホームオスタティックな睡眠需要（sleep pressure）」と関係する指標として解釈されます。一方、覚醒時に前頭部デルタが増える場合、病的な変化や機能低下のサインとして議論されることもあります。

  • シータ（4–8 Hz 前後） ヒトの頭皮 EEG では、眠気・うとうとした状態で増えることが多く、覚醒度の低下の指標として使われます。 認知神経科学では、前頭部シータ活動がワーキングメモリや認知制御、エラー検出、学習などと関係することが多数報告されており、「単なる眠気の指標」ではなく、認知負荷や記憶処理と結びつくリズムとしても解釈されています。

  • アルファ（8–12/13 Hz 前後） 後頭部を中心に、閉眼安静状態で最もよく観察される律動で、いわゆる「後頭部アルファ」が典型例です。 伝統的には「安静・リラックス」の指標と説明されることが多いですが、近年の研究では、アルファ活動はむしろ感覚入力や情報処理を抑制する“トップダウン抑制”のメカニズムと関連づけて解釈されることが増えています。 たとえば、視覚課題で注意を向けていない側の視野に対応する頭皮部位ではアルファパワーが増加し、「その領域の処理が抑えられている」ことを反映する、といった結果が多く報告されています。

  • ベータ（13–30 Hz 前後） 一般には、覚醒・集中・運動準備などと関係づけられます。運動野周辺では、安静時にベータ活動が強く、実際に運動するときや運動のイメージ時にベータパワーが低下する（event-related desynchronization, ERD）ことが知られています。 臨床的には、パーキンソン病などの運動障害における異常なベータ同期が議論されており、深部脳刺激（DBS）のターゲット選択とも関係しています。

  • ガンマ（30 Hz 以上） 高周波数の活動で、注意、知覚の統合、作業記憶、意識状態などと関連づけた研究が多数あります。特定の刺激や課題条件で局所的にガンマ活動が増えることが報告されており、「局所回路の同期」や「情報統合」の指標として解釈されることが多いです。 一方で、頭皮 EEG の高周波成分には筋電図（EMG）による汚染が混ざりやすいという問題もあり、ガンマ活動の解釈には注意が必要です。最近の研究では、脳表（ECoG）や深部記録、MEG など、筋電の影響をより低減できる手法での検証が進んでいます。

■ 2.3 分野による定義の違いと最新の動向

　近年の研究動向としては、次のような傾向があります。

  • 睡眠研究 睡眠段階ごとの特徴を扱うため、デルタ帯を 0.5–4 Hz よりも狭く定義したり、0.5–1 Hz の「超低周波成分（slow oscillation）」を特に区別して扱う研究もあります。また、年齢や発達段階によってアルファ帯のピーク周波数（individual alpha frequency, IAF）が変化することから、個人ごとのピーク周波数に合わせてバンド定義を調整する試みも行われています。

  • 認知神経科学・脳機能イメージングと組み合わせた EEG/MEG アルファ・シータ・ベータなどを一律に固定周波数帯で扱うのではなく、

    • 個人ごとのピーク周波数
    • リズムごとの空間分布（どの脳領域で強いか） を考慮して、「機能単位としてのリズム」をより柔軟に定義しようとする流れがあります。また、ガンマ帯は 30–45 Hz 程度までの「低ガンマ（low gamma）」と、70–100 Hz 近辺の「高ガンマ（high gamma）」とに分けて扱う研究が増えています。

  • 臨床 EEG 臨床分野では、周波数バンドの解釈に加えて、波形の形（スパイク・鋭波・徐波など）や、空間的な分布、発生パターンが重視されます。バンドパワーそのものよりも、「どの部位でどのような異常波が出ているか」という視点が優先される場合も多く、研究用途とは少し違った言葉づかいをすることがあります。

■ 2.4 ここで採用するバンド定義

　以上のように、周波数バンドの定義は分野によって多少異なりますが、本記事では ヒト頭皮 EEG の解析で広く用いられている、比較的無難な定義として、

  • Delta：0.5–4 Hz
  • Theta：4–8 Hz
  • Alpha：8–13 Hz
  • Beta：13–30 Hz
  • Gamma：30–45 Hz

を採用します。この帯域は、睡眠研究・覚醒時の認知課題・基本的な脳波解析において、まず押さえておくべき「オーソドックスな設定」のひとつといえます。

3. Rでやってみる

　脳波解析では、長時間の脳波をまとめて扱うのではなく、一定の長さごとに区切って解析することがよく行われます。この区切ったひとまとまりを エポック（epoch） と呼びます。たとえば「30 秒エポック」といえば、記録された脳波を 30 秒ずつに分けましたということです。
　エポックはとくに睡眠研究（睡眠段階判定）や長時間 EEG の解析で標準的な区切り方ですが、認知課題や臨床 EEG でも、「時間経過による変化」を追うために広く用いられています。 　ここでは、1 エポックごとに Welch 法でパワースペクトルを推定し、デルタ・シータ・アルファ・ベータ・ガンマといった脳波バンドのパワーを計算します。これにより、脳の状態が時間とともにどのように変化したかを、定量的に評価できます。

■ 今回のスクリプトの概要

　ここから示す R スクリプトでは、EDF から読み取った脳波信号 ch1 を使って、

  1. 30 秒エポックに分割し
  2. 各エポックで Welch 法により PSD（パワースペクトル密度）を計算し
  3. 脳波バンドごとのパワー（デルタ・シータ・アルファ・ベータ・ガンマ）を求め
  4. エポックの時間情報（絶対時刻）と一緒にテーブル化

するところまで自動的に行います。
　最初にも述べた通り、今回の前提として扱う EDF の要素は以下の 3 つです。

sig <- ch1$signal        # EEG 波形ベクトル
fs  <- ch1$sRate         # サンプリング周波数 [Hz]
start_time0 <- as.POSIXct(ch1$startTime)  # 計測開始時刻

これらの情報がそろっていれば、長時間脳波を「30 秒ごとの脳波指標」に変換することができます。 以下に、実際の R スクリプトを示します。

############################################################
# EDF から読み込んだ EEG (ch1) について，
# 30 秒エポックごとに Welch 法で PSD を推定し，
# 各エポックの脳波バンドパワー（Δ・θ・α・β・γ）を求めるスクリプト
#
# 前提：
#   ch1$signal    : EEG 波形ベクトル（単位は μV など）
#   ch1$sRate     : サンプリング周波数 [Hz]
#   ch1$startTime : 記録開始時刻（POSIXlt / POSIXct / 文字列 いずれでも可）
#
# 出力：
#   EEG_band_power : data.frame
#     - epoch          : 30 秒エポックの通し番号（1, 2, 3, ...）
#     - center_sec     : 記録開始からの経過時間（エポック中心）[s]
#     - center_time    : 実際の時計時刻（エポック中心，startTime がある場合）
#     - P_delta, ...   : 各脳波バンドのパワー（絶対値，単位^2）
#     - P_total_0.5_45Hz : 0.5–45 Hz の全パワー（絶対値）
############################################################


## 1. 必要な情報を取り出す ----------------------------------------

# EEG 波形（1 チャンネル分）
sig <- ch1$signal          # 例：μV 単位の連続データ

# サンプリング周波数 [Hz]
fs  <- ch1$sRate           # 例：250, 256, 512, 1000 など

# 記録開始時刻（なければ NA にしておく）
# ch1$startTime が文字列でも POSIX* でも動くように tryCatch でラップ
start_time0 <- tryCatch(
  as.POSIXct(ch1$startTime),
  error = function(e) NA
)

# 全サンプル数（＝記録された点の数）
N <- length(sig)

# 30 秒エポックの長さ（サンプル数）
epoch_sec <- 30                               # エポックの長さ [s]
epoch_len <- as.integer(round(epoch_sec * fs))# 30 秒に相当するサンプル数

# 取れるエポック数（端数は切り捨て）
n_epoch <- floor(N / epoch_len)
if (n_epoch < 1) {
  stop("30秒エポックが 1 つも取れません。記録時間が短すぎます。")
}


## 2. Welch 法の設定 -----------------------------------------------

# Welch 法で用いる短い窓（セグメント）の長さを決める
# ここでは「5 秒窓，50% オーバーラップ」を例として用いる
win_sec  <- 5
win_len  <- as.integer(round(win_sec * fs))   # 5 秒分のサンプル長
step_len <- as.integer(floor(win_len / 2))    # 50% オーバーラップ → 窓の半分ずつずらす

# 窓関数（ここでは Hann 窓 / Hanning window）
# 端の不連続性をなめらかにして，スペクトル漏れ（leakage）を軽減する目的
hann_window <- 0.5 - 0.5 * cos(2 * pi * (0:(win_len - 1)) / (win_len - 1))


## 3. Welch PSD を計算する関数 -------------------------------------
# 与えられた 1 エポック分の波形 x について，
#   - Hann 窓をかけた短いセグメントに分割し
#   - 各セグメントごとに FFT をとってペリオドグラムを計算し
#   - それらを平均して Welch PSD を推定する
# 返り値は 0 ～ ナイキスト周波数までの
#   list(freq = 周波数ベクトル, psd = PSD ベクトル)

welch_psd <- function(x, fs, win_len, step_len, window) {
  # x        : 1 エポック分の波形（数値ベクトル）
  # fs       : サンプリング周波数 [Hz]
  # win_len  : 窓長（サンプル数）
  # step_len : セグメントのシフト長（サンプル数）
  # window   : 窓関数ベクトル（長さ win_len）

  n <- length(x)
  if (n < win_len) {
    stop("エポック長が窓長より短いです。win_sec を短くするか，epoch_sec を長くしてください。")
  }

  # 窓の開始インデックス列（1 始まり）
  # 例：win_len = 500, step_len = 250 のとき
  # 1, 251, 501, ... のようにずらしていく
  starts <- seq(1, n - win_len + 1, by = step_len)
  n_win  <- length(starts)        # 使用されるセグメント数

  # FFT の周波数軸（0 ～ (fs * (win_len - 1) / win_len)）
  freq <- (0:(win_len - 1)) * fs / win_len

  # 片側スペクトル用のインデックス（0 ～ ナイキストまで）
  # たとえば win_len が偶数なら 0, 1, ..., fs/2 まで
  idx_half <- 1:(floor(win_len / 2) + 1)

  # PSD を加算していく器（0 で初期化）
  psd_accum <- numeric(length(idx_half))

  # 窓関数のエネルギー（|window|^2 の総和）
  # 正しい PSD (単位^2 / Hz) に正規化するために使う
  win_norm <- sum(window^2)

  # --- メインループ：各セグメントごとにペリオドグラムを加算 ---
  for (k in seq_len(n_win)) {

    # このセグメントに対応する x のインデックス
    idx <- starts[k]:(starts[k] + win_len - 1)

    # セグメント波形を取り出す
    seg <- x[idx]

    # DC オフセット（平均値）を除去
    # 低周波成分に偏りがあるとスペクトル推定が歪むため
    seg <- seg - mean(seg)

    # 窓関数をかける（端のなまりをつける）
    segw <- seg * window

    # FFT（R の fft() は 0～N-1 の周波数成分を返す）
    # 正規化は後で手動で行う
    Xk <- fft(segw)

    # 二側スペクトルのパワー |X(f)|^2 を計算し，
    # サンプリング周波数と窓のエネルギーで割ることで PSD に正規化する
    # 単位は「(元信号の単位)^2 / Hz」
    Pk <- (Mod(Xk)^2) / (fs * win_norm)

    # 0 ～ ナイキスト周波数の部分だけを取り出して加算
    psd_accum <- psd_accum + Re(Pk[idx_half])
  }

  # セグメントごとに平均（Welch 法の「平均をとる」部分）
  psd <- psd_accum / n_win

  # 片側 PSD にするため，
  # DC 成分（0 Hz）とナイキスト周波数を除く部分は 2 倍する
  # （元の二側スペクトルの両側を片側に集約するイメージ）
  if (length(psd) > 2) {
    psd[2:(length(psd) - 1)] <- 2 * psd[2:(length(psd) - 1)]
  }

  # 0 ～ ナイキスト周波数の周波数軸
  freq_half <- freq[idx_half]

  # 結果を返す
  list(freq = freq_half, psd = psd)
}


## 4. バンドパワーを計算する関数 -----------------------------------
# 与えられた PSD (freq, psd) から，
#   [f_low, f_high) の周波数帯に含まれるパワー（積分値）を求める
# 離散周波数なので，和 × 周波数分解能 df で近似する

band_power <- function(freq, psd, f_low, f_high) {
  # freq : 周波数軸（0 ～ ナイキスト）
  # psd  : PSD 値（freq と同じ長さ）
  # f_low, f_high : 積分する周波数帯 [Hz]

  # 対象帯域に含まれるインデックス
  idx <- which(freq >= f_low & freq < f_high)
  if (length(idx) == 0) {
    return(NA_real_)   # 対象帯域が存在しない場合は NA
  }

  # 周波数分解能（隣り合う周波数ビンの間隔）
  # ※等間隔を前提としている
  df <- freq[2] - freq[1]

  # PSD の積分 ≒ PSD の和 × df
  sum(psd[idx]) * df
}


## 5. 脳波バンドの定義 ---------------------------------------------
# 解析目的に応じて境界周波数は調整してください。
# ここでは一例として：
#   delta : 1–4 Hz
#   theta : 4–8 Hz
#   alpha : 8–13 Hz
#   beta  : 13–30 Hz
#   gamma : 30–45 Hz

bands <- list(
  delta = c(0.5, 4),
  theta = c(4, 8),
  alpha = c(8, 13),
  beta  = c(13, 30),
  gamma = c(30, 45)
)

band_names <- names(bands)


## 6. 各 30 秒エポックで PSD & バンドパワーを計算 -------------------

# 結果を格納する器をあらかじめ用意しておく

# 各エポックの中心時刻（記録開始からの経過時間 [s]）
epoch_center_sec  <- numeric(n_epoch)

# 各エポックの中心時刻（実際のカレンダー時刻）
epoch_center_time <- rep(as.POSIXct(NA), n_epoch)

# 各エポック × 各バンドのパワーを格納する行列
band_mat <- matrix(NA_real_,
                   nrow = n_epoch,
                   ncol = length(bands))

# 各エポックの「全パワー」（ここでは 0.5–45 Hz の帯域とする）
total_pow <- numeric(n_epoch)

# --- エポックごとにループ ---
for (e in seq_len(n_epoch)) {

  # このエポックに対応するサンプルの開始・終了インデックス
  idx_start <- (e - 1) * epoch_len + 1
  idx_end   <- idx_start + epoch_len - 1

  # 1 エポック分の EEG 波形を切り出す
  x_epoch <- sig[idx_start:idx_end]

  # このエポックについて Welch PSD を推定
  wpsd <- welch_psd(x_epoch, fs, win_len, step_len, hann_window)
  freq <- wpsd$freq
  psd  <- wpsd$psd

  # 各脳波バンドごとにバンドパワーを計算
  for (b in seq_along(bands)) {
    f_band <- bands[[b]]   # 例：c(4, 8) など
    band_mat[e, b] <- band_power(freq, psd, f_band[1], f_band[2])
  }

  # 全パワー（例として 0.5–45 Hz）
  total_pow[e] <- band_power(freq, psd, 0.5, 45)

  # このエポックの中心時刻（記録開始からの経過秒）
  # e = 1 → 15 秒, e = 2 → 45 秒, ... というイメージ
  center_sec <- (e - 0.5) * epoch_sec
  epoch_center_sec[e] <- center_sec

  # 記録開始時刻がわかっていれば，絶対時刻も求める
  if (!is.na(start_time0)) {
    epoch_center_time[e] <- start_time0 + center_sec
  }
}


## 7. 結果テーブルを作成 --------------------------------------------

# 基本情報（エポック番号と時刻）
EEG_band_power <- data.frame(
  epoch       = 1:n_epoch,
  center_sec  = epoch_center_sec,
  center_time = epoch_center_time
)

# バンドパワーの列を追加
colnames(band_mat) <- paste0("P_", band_names)  # 例：P_delta, P_theta, ...
EEG_band_power <- cbind(EEG_band_power, band_mat)

# 全パワー列を追加（ここでは 0.5–45 Hz）
EEG_band_power$P_total_0.5_45Hz <- total_pow


## 8. 結果を表示／保存 ----------------------------------------------

# 結果の先頭数行を確認
head(EEG_band_power)

# 必要であれば CSV に保存する
# write.csv(EEG_band_power,
#           "EEG_band_power_30sWelch.csv",
#           row.names = FALSE)

　このスクリプトは「30 秒エポック × Welch 法 × バンドパワー抽出」という、一連の EEG 時系列解析の基本パイプラインをそのまま実装しています。ここでは、コードの中で特に重要な部分を分解しながら説明します。

■ エポック分割（30 秒）

　エポックの長さ（秒単位）は以下の部分で指定しています。

epoch_sec <- 30
epoch_len <- as.integer(round(epoch_sec * fs))

睡眠解析や長時間モニタリングの世界では 30 秒エポックが標準的によく使われています。ただし、研究目的によっては、

  • 1 秒〜5 秒エポック（認知課題）
  • 60 秒エポック（粗い時間変化を見る）
  • イベントに合わせた可変長エポック

などに変えることもあります。その場合は単に epoch_sec を変更すればよい構造になっています。

■ Welch 法の設定（窓長 5 秒 + 50% オーバーラップ）

　Welch 法で用いる部分時系列の長さ（秒単位）は、以下の部分で指定しています。

win_sec  <- 5
step_len <- floor(win_len / 2)

Welch 法は、「窓長（win_sec）」、「オーバーラップ（ここでは 50%）」の設定に強く依存します。
　窓長を短くすれば時間分解能が上がり、窓長を長くすれば周波数分解能が上がります。
一般的には

  • 1～8 秒窓がよく使われる
  • 50% オーバーラップが無難（慣例的）
  • 心電図・脳波のような生体信号では 4～8秒窓が多い

という傾向があります。

■ Hann 窓（ハニング窓）の使用

　部分時系列の切り出しには、ハニング窓を使っています。その定義は以下の部分です。

hann_window <- 0.5 - 0.5 * cos(...)

Hann 窓はスペクトル漏れ（leakage）を抑える効果が高く、Welch 法で使われることが多い窓です。他の窓関数を使っても構いませんが、脳波解析では Hann は無難な選択です。

■ 正しく PSD を計算するポイント

　Welch 法の実装では、いくつか重要な正規化があります。

  • FFT の結果をそのまま使わず、窓のエネルギー sum(window^2）で割る
  • 周波数 0〜ナイキストまでの片側 PSDを作るため、 0 Hz とナイキスト周波数以外を 2 倍する
  • 周波数分解能 df を計算し、PSD の積分値（バンドパワー）を求めるときに df を掛ける

これらは、「物理単位を正しく保った PSD を得るために必須」の処理です。

■ バンド定義は研究目的で調整する

　以下の部分で、バンドの周波数帯 (Hz)を定義しています。

bands <- list(
  delta = c(0.5, 4),
  theta = c(4, 8),
  alpha = c(8, 13),
  beta  = c(13, 30),
  gamma = c(30, 45)
)

バンドの定義は絶対ではありません。分野・研究室・目的によって異なります。

4. まとめ

　大事なのは、まずは実際にスクリプトを動かしてみることです。自分のデータで走らせてみて、「この被験者はデルタが高いな」「覚醒するとアルファが減ってベータが増えているな」といった感覚を、自分の目で確かめてみてください。とにかく手を動かして試してみること。そして、その結果から「なぜそうなるのか」を考えていくこと。そうした試行錯誤の時間こそが、脳波解析の感覚を自分のものにしていく、いちばん確かなプロセスだと思います。　　

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。

ランキング参加中
知識
ランキング参加中
テクノロジー
ランキング参加中
数学・科学・工学