Rで解析：頭蓋内脳波の信号処理と3次元データをプロット！「ravetools」パッケージの紹介

頭蓋内電極を用いた脳波（iEEG）の記録では、数分から数時間に及ぶ高解像度の信号を扱うことになります。しかし、商用電源に由来するノイズの除去や時間周波数解析を長時間の信号に対して行うには、計算時間とメモリの両面で手間がかかります。

「ravetools」パッケージは、iEEG信号の処理を目的とした関数群を収録しています。ノッチフィルタやバンドパスフィルタの設計と適用、FFTW3を利用した高速フーリエ変換、Welchペリオドグラムやウェーブレット変換による時間周波数解析が可能です。また、3次元からの等値面の生成やメッシュの描画、表面上の最短距離の計算といった空間データの処理のコマンドも収録されています。本パッケージの利用で、長時間の信号処理と3次元データの確認を同一の環境で進められるのではないかと考えます。

パッケージバージョンは0.3.0。Windows 11 x64 (build 26200)のR version 4.6.1で確認しています。

＜おすすめのRに関する書籍です＞

Rによる統計データ解析 | 小池 祐太, 村田 昇, 吉田 朋広

確率論の基礎から推定・検定・回帰までを、数式とRのコードの両方で追える大学教科書。理論をきちんと押さえたい人向け。

パッケージのインストール

下記コマンドを実行してください。

# パッケージのインストール
install.packages("ravetools")

# パッケージの読み込み
library("ravetools")

コマンド例

詳細はコメント、パッケージのヘルプを確認してください。

iEEGは、頭蓋内に留置した電極から記録する脳波を指します。FFTは高速フーリエ変換（Fast Fourier Transform）の略称で、時間の経過に沿って並ぶ信号を周波数の成分へ分解する計算手法です。

以降のコマンド例では、架空の観測信号と、北海道の地点名を用いた架空の観測点網を使用します。はじめに、記事全体で使用する信号データを作成します。

# 乱数の種を固定
set.seed(20260701)

# サンプリング周波数（Hz）を設定
sample_rate <- 1000

# 記録時間（秒）を設定
duration <- 10
# 時間軸を作成
time_sec <- seq(0, duration, by = 1 / sample_rate)

# 観測信号を作成（8Hzの徐波、40Hzの速波、50Hzの商用電源ノイズとその第2高調波、白色雑音）
signal_raw <- 30 * sin(2 * pi * 8 * time_sec) +
  12 * sin(2 * pi * 40 * time_sec) +
  20 * sin(2 * pi * 50 * time_sec) +
  8 * sin(2 * pi * 100 * time_sec) +
  rnorm(length(time_sec), sd = 5)

デジタルフィルタの診断：diagnose_filterコマンド

設計したデジタルフィルタの周波数応答と、サンプル信号へ適用した際の挙動を並べて描画します。通過帯域が意図どおりかを確認する際に使用します。

# ナイキスト周波数を計算
nyquist <- sample_rate / 2

# 通過帯域（1〜40Hz）を正規化周波数へ変換
pass_band <- c(1, 40) / nyquist

# FIRフィルタを次数200で設計（移動平均型のため自己回帰係数は1）
fir_filter <- fir1(200, pass_band, "pass")

# 設計したFIRフィルタの周波数応答を診断
diagnose_filter(
  b = fir_filter$b, a = fir_filter$a, fs = sample_rate,
  n = 1024, vlines = c(1, 40), cutoffs = c(-3, -6, -12)
)

# バターワースフィルタを次数4で設計
iir_filter <- butter(4, pass_band, "pass")

# 設計したバターワースフィルタの周波数応答を診断
diagnose_filter(
  b = iir_filter$b, a = iir_filter$a, fs = sample_rate,
  n = 1024, vlines = c(1, 40), xlim = "auto"
)

位相の遅れを伴わないフィルタ処理：filtfiltコマンド

信号に対して順方向と逆方向の両側からフィルタを適用します。片方向のみの処理で生じる位相の遅れが相殺されるため、波形の時間的な位置を保ったまま帯域を制限できます。

# 1〜40Hzを通過帯域とするバターワースフィルタを設計
band_filter <- butter(4, c(1, 40) / (sample_rate / 2), "pass")

# 順方向と逆方向の両側からフィルタを適用
signal_zerophase <- filtfilt(band_filter$b, band_filter$a, signal_raw)

# 出力の長さが入力と一致することを確認
length(signal_zerophase) == length(signal_raw)
[1] TRUE

表面上の最短距離の計算：dijkstras_surface_distanceコマンド

ダイクストラ法により、起点となるノードから全ノードまでの最短距離を計算します。facesが2列の場合はグラフの辺、3列の場合は3次元メッシュの三角形として扱われます。

# 観測点8地点の平面座標（東西方向と南北方向の距離、単位はキロメートル）を作成
station_position <- matrix(
  c(  0,   0,
      8,  -9,
     -6,   9,
    -18,  22,
     14, -26,
     26,  14,
      6,  34,
     52,  46),
  ncol = 2, byrow = TRUE)

# 各観測点の名称を設定
station_name <- c("恵庭", "千歳", "北広島", "札幌",
                  "苫小牧", "夕張", "岩見沢", "富良野")

# 観測点どうしを結ぶ経路を2列の行列として作成
station_edge <- matrix(
  c(1, 2,
    1, 3,
    3, 4,
    2, 5,
    2, 6,
    4, 7,
    3, 7,
    6, 7,
    6, 8,
    7, 8),
  ncol = 2, byrow = TRUE)

# 恵庭（1番）を起点として全観測点までの最短距離を計算
station_distance <- dijkstras_surface_distance(
  positions = station_position,
  faces = station_edge,
  start_node = 1,
  face_index_start = 1)

# 計算結果に含まれる要素を確認
names(station_distance)
[1] "paths"               "start_node"          "face_index_start"
[4] "max_search_distance" "n_nodes"             "n_faces"

# 対象となったノード数と経路数を確認c(station_distance$n_nodes, station_distance$n_faces)
[1]  8 10

＜おすすめのRに関する書籍です＞

Rで楽しむ統計 (Wonderful R 1)

検定や信頼区間、回帰といった基本を、Rで手を動かしながら確認していく一冊。

www.amazon.co.jp

最短経路の取得：surface_pathコマンド

dijkstras_surface_distance()の計算結果から、起点と目標ノードを結ぶ最短経路を取り出します。経路上のノード番号と、起点からの累積距離を持つデータフレームが返されます。

# 恵庭から富良野（8番）までの最短経路を取得
station_path <- surface_path(station_distance, target_node = 8)

# 経路上のノード番号と起点からの累積距離を表示
station_path
  path distance
1    1  0.00000
2    2 12.04159
3    6 41.24776
4    8 82.47881

# ノード番号を観測点の名称へ置き換えて表示
station_name[station_path$path]
[1] "恵庭"   "千歳"   "夕張"   "富良野"

実数信号の高速フーリエ変換：fftw_r2cコマンド

実数の信号を複素数のスペクトルへ変換します。HermConjに1を指定すると入力と同じ長さの完全なスペクトルが、0を指定すると重複のない片側スペクトルが返されます。

# 解析対象として先頭の1024点を切り出し
signal_head <- signal_raw[seq_len(1024)]

# 完全なエルミート対称スペクトルを取得
spectrum_full <- fftw_r2c(signal_head, HermConj = 1)

# 標準のfftによる結果と一致することを確認
all.equal(spectrum_full, stats::fft(signal_head))
[1] TRUE

# 重複のない片側スペクトルを取得
spectrum_half <- fftw_r2c(signal_head, HermConj = 0)
# 片側スペクトルの長さを確認
length(spectrum_half)
[1] 513

# 直流成分を除いて振幅が最大となる位置を求める
peak_index <- which.max(Mod(spectrum_half[-1])) + 1

# 位置を周波数（Hz）へ換算
(peak_index - 1) * sample_rate / length(signal_head)
[1] 7.8125

1024点を1000Hzで切り出しているため周波数の分解能は約0.977Hzとなり、8Hzの成分は7.8125Hzの位置に現れます。

複素数信号の高速フーリエ変換：fftw_c2cコマンド

複素数の入力に対して順方向と逆方向の変換を行います。逆方向の変換は正規化されないため、元の信号へ戻す際は要素数で割る必要があります。

# 実部と虚部を持つ複素数の信号を作成
signal_complex <- complex(real = signal_head, imaginary = rev(signal_head))

# 順方向の変換を実行
spectrum_c2c <- fftw_c2c(signal_complex, inverse = 0)

# 標準のfftによる結果と一致することを確認
all.equal(spectrum_c2c, stats::fft(signal_complex))
[1] TRUE

# 逆方向の変換を実行し要素数で割って正規化
signal_restored <- fftw_c2c(spectrum_c2c, inverse = 1) / length(signal_complex)

# 復元した信号が元の信号と一致することを確認
all.equal(signal_restored, signal_complex)
[1] TRUE

列ごとの一括変換：mvfftw_r2cコマンド

行列の各列を独立して変換します。多チャネルの信号をまとめて処理する際に使用します。fftw_r2c()とは引数の並び順が異なり、HermConjが3番目に置かれ、初期値も0Lである点に注意が必要です。

# 4チャネル分の信号を列方向に並べた行列を作成
channel_matrix <- matrix(
  c(signal_head,
    signal_head * 0.8 + rnorm(1024, sd = 3),
    signal_head * 0.6 + rnorm(1024, sd = 3),
    signal_head * 0.4 + rnorm(1024, sd = 3)),
  ncol = 4
)

# 完全なスペクトルを指定して列ごとに変換
spectrum_matrix <- mvfftw_r2c(channel_matrix, fftwplanopt = 0, HermConj = 1)
# 標準のmvfftによる結果と一致することを確認
all.equal(spectrum_matrix, stats::mvfft(channel_matrix + 0i))
[1] TRUE

# 出力の次元を確認
dim(spectrum_matrix)
[1] 1024    4

＜おすすめのRに関する書籍です＞

はじめてのR: ごく初歩の操作から統計解析の導入まで | 村井 潤一郎

本当に最初の一歩から始めるR入門。インストールと基本操作から統計解析の入口までを扱う。

メッシュの点群によるプロット：plot_mesh_dotcloudコマンド

メッシュの頂点を点の集合として描画します。形状の密度や分布を把握したい場合や、面を持たない電極の座標を重ねて表示したい場合に適しています。

# 3次元ボリュームの一辺の分割数を設定
grid_n <- 48

# 各軸の座標値を作成
axis_seq <- seq(-1.5, 1.5, length.out = grid_n)

# x座標を3次元配列として展開
coord_x <- array(axis_seq, dim = c(grid_n, grid_n, grid_n))

# 軸を入れ替えてy座標を作成
coord_y <- aperm(coord_x, c(2, 1, 3))

# 軸を入れ替えてz座標を作成
coord_z <- aperm(coord_x, c(3, 2, 1))

# 円環（トーラス）形状を表す距離場を計算
torus_field <- (sqrt(coord_x^2 + coord_y^2) - 0.9)^2 + coord_z^2

# 距離場を二値のボリュームデータへ変換
torus_volume <- array(as.numeric(torus_field < 0.25^2), dim = dim(torus_field))

# ボリュームデータから等値面のメッシュを生成
torus_mesh <- vcg_isosurface(torus_volume, threshold_lb = 0.5)

# 生成したメッシュを点群としてプロット
plot_mesh_dotcloud(
  torus_mesh,
  eye = c(0, -120, 90),
  lookat = c(0, 0, 0),
  up = c(0, 0, 1),
  col = "steelblue",
  cex = 1.5)

# メッシュの重心を求める
mesh_centre <- rowMeans(torus_mesh$vb[1:3, ])

# 電極を模した点の数を設定
n_electrode <- 24

# 円周上へ等間隔に配置するための角度を作成
electrode_angle <- seq(0, 2 * pi, length.out = n_electrode + 1)[-(n_electrode + 1)]

# 電極の座標をmesh3d形式のリストとして構成
electrode_mesh <- structure(
  list(vb = rbind(
    mesh_centre[1] + 14 * cos(electrode_angle),
    mesh_centre[2] + 14 * sin(electrode_angle),
    mesh_centre[3] + rnorm(n_electrode, sd = 1.5)
  )),
  class = "mesh3d")

# 表面のメッシュと電極の点群を同時にプロット
plot_mesh_dotcloud(
  mesh = list(torus_mesh, electrode_mesh),
  eye = c(0, -120, 90),
  lookat = c(0, 0, 0),
  up = c(0, 0, 1),
  col = list("steelblue", "orangered"),
  pch = c(16L, 17L),
  cex = c(1.5, 2.5))

・生成したメッシュを点群としてプロット

・表面のメッシュと電極の点群を同時にプロット

メッシュのポリゴンのプロット：plot_mesh_polygonコマンド

メッシュの三角形をランバートシェーディングで塗り分けて描画します。表面の凹凸を陰影として確認したい場合に適しています。面を持たないメッシュは、頂点ごとに半径cexの球へ置き換えてプロットします。

### 注意：plot_mesh_dotcloudコマンドの実行例を実施してから実行してください。
# 円環のメッシュをポリゴンとしてプロット
plot_mesh_polygon(
  torus_mesh,
  eye = c(0, -120, 90),
  lookat = c(0, 0, 0),
  up = c(0, 0, 1),
  col = "steelblue",
  main = "円環メッシュの陰影表示")

# 表面を半透明にし、電極を球へ置き換えて重ねてプロット
plot_mesh_polygon(
  mesh = list(torus_mesh, electrode_mesh),
  eye = c(0, -120, 90),
  lookat = c(0, 0, 0),
  up = c(0, 0, 1),
  col = list("steelblue", "orangered"),
  alpha = c(0.45, 1),
  cex = 1.8,
  sphere_subdivision = 2L,
  ambient_intensity = 0.3)

・円環のメッシュをポリゴンとしてプロット

・表面を半透明にし、電極を球へ置き換えて重ねてプロット

チャネル信号の診断プロット：diagnose_channelコマンド

1本または2本の信号について、波形、Welchペリオドグラム、ヒストグラムを並べた診断用の図を作成します。フィルタ処理の前後を比較する際に使用します。

# 北海道の商用電源周波数50Hzとその高調波を対象にノッチフィルタを適用
signal_notched <- notch_filter(
  signal_raw, sample_rate,
  lb = c(49, 98, 148),
  ub = c(51, 102, 152)
)

# 処理の前後を並べて診断用の図を作成
diagnose_channel(
  signal_raw, signal_notched,
  srate = sample_rate,
  name = c("原信号", "50Hz除去後"),
  max_freq = 200,
  cex = 1,
  main = "恵庭観測点の信号診断"
)

複数信号の同時プロット：plot_signalsコマンド

行列の各行を1本のトレースとして、縦方向に間隔を空けて並べて描画します。多チャネルの信号を時間軸に沿って比較する際に使用します。

# チャネル数を設定
n_channel <- 4
# チャネルごとに位相をずらした信号を行方向に並べた行列を作成
electrode_signals <- t(vapply(seq_len(n_channel), function(i) {
  30 * sin(2 * pi * 8 * time_sec + i * pi / 3) +
    10 * sin(2 * pi * 40 * time_sec) +
    rnorm(length(time_sec), sd = 4)
}, FUN.VALUE = numeric(length(time_sec))))
# 各チャネルの名称を設定
electrode_names <- c("恵庭-01", "恵庭-02", "千歳-01", "千歳-02")

# 全チャネルを縦に並べて描画
plot_signals(
  electrode_signals,
  sample_rate = sample_rate,
  channel_names = electrode_names
)

# 開始から2秒経過した時点より3秒分のみを切り出して描画
plot_signals(
  electrode_signals,
  sample_rate = sample_rate,
  channel_names = electrode_names,
  start_time = 2,
  duration = 3,
  col = c("steelblue", "orangered")
)

＜おすすめのRに関する書籍です＞

Rで楽しむベイズ統計入門［しくみから理解するベイズ推定の基礎］ Data Science Library

ベイズ推定のしくみを、仕組みの理解から丁寧に追う入門書。

www.amazon.co.jp

この記事が誰かの役に立ちますように。