理論から実装まで
Morlet、CWT/DWT、COI、phase、cluster permutation、PyWavelets、歩行EMG・運動学を切れ目なく扱います。
REVIEW & REDESIGN
0. まず、教材をどう読みやすくしたか
元資料は内容の幅が広く、実装に必要な注意点も多く含んでいました。今回は「知識を増やす」より「迷子にならない順番に置く」ことを優先しました。
問いを先に置く
数式やコードを先に暗記せず、「何を測りたいか」と「何は言えないか」を各章の先頭に置きました。
歩行研究の判断表
実時間とgait-cycle正規化、strideと被験者、powerと生理学的解釈を分けて考えられるようにしました。
このページを読み終えたら、説明できることと、解析コードを書き始める前に決めることが分かっている。
LEARNING MAP
1. 6つの階段で理解する
各ステップで「一文で説明できる」状態を作ってから、次へ進みます。
1時間–周波数なぜFourierだけでは足りない？2局所waveletGaussianと振動を重ねる3係数W内積＝似ている度合い4power / phase強さと周期上の位置5地図を読む帯・島・ridge・COI6研究へ接続問い・単位・統計・解釈
学習のコツ：「scaleは周波数そのもの」ではありません。waveletの伸び縮みを表すパラメータで、実装では対応する周波数へ変換して読む、と分けて覚えます。
CHAPTER 1 · WHY
2. waveletは「いつ」を失わない
歩行のように、同じ記録の中で短いイベントと周期的な動きが共存する信号では、周波数だけでなく時間位置が重要です。
| 方法 | 得意な問い | 窓の考え方 | 歩行研究での使いどころ |
|---|---|---|---|
| Fourier | 記録全体に何Hzがある？ | 時間全体をまとめて見る | 全試行の主な周波数を把握 |
| STFT | 固定幅の窓で、いつ何Hz？ | すべての周波数で窓幅が同じ | 条件をそろえた局所解析 |
| Wavelet / CWT | どの時間スケールの成分がいつ？ | 高周波は短く、低周波は長く見る | burst、impact、歩行周期の共存 |
最低限の式
x(t) → W(a,b) = ∫ x(t) ψ*a,b(t) dt
a はscale（伸び縮み）、b はwaveletの中心を置く時間位置、* は複素共役です。
時間–周波数のトレードオフ
短い窓：「いつ」は細かいが、近い周波数を区別しにくい。
長い窓：周波数は細かいが、イベント時刻がぼやける。
waveletはこの窓幅を周波数に応じて変える設計です。
Fourier＝何Hz、wavelet＝何Hzがいつ。歩行の「いつ」を扱う理由を一文で言える。
CHAPTER 2 · MORLET
3. Morlet waveletを触ってみる
Morletは、振動する波をGaussianの封筒で包んだ「局所的なテンプレート」です。下のスライダーでscaleと時間位置を動かします。
インタラクティブ：waveletの伸び縮みと移動
上段は信号、下段は選択中のwavelet。青い山がGaussian envelope、紫の波が局所振動です。
入力信号 x(t)局所振動Gaussian envelope
scale a
a ↑ → waveletが広がる → 対応周波数は低くなる。
a ↓ → waveletが縮む → 対応周波数は高くなる。
周期とGaussianの幅が一緒に変わるので、周波数と時間分解能を別々に操作しているわけではありません。
中心 b
b は「位相」ではなく、waveletの中心を置く時間位置です。
理論上は連続値ですが、デジタル実装では時間点とscale候補を細かく走査して近似します。
scale a＝伸び縮み、b＝中心位置。低周波用waveletほど時間幅が長い、と図を見ながら説明できる。
CHAPTER 3 · COEFFICIENT
4. W(a,b)は「その形がどれだけあるか」
各scaleと各時間位置にwaveletを置き、入力信号との重なりを計算します。形とタイミングが合えば、内積の大きさが大きくなります。
① waveletを選ぶscale a と中心 b を決める
② 掛けて足すx(t) と ψa,b(t) の内積
③ 係数を得るW(a,b) が大きいほど似ている
PCAとの類推
PCA: score = x · w
Wavelet: W(a,b) = x · ψa,b
Wavelet: W(a,b) = x · ψa,b
どちらも「元データを、あるパターンへ射影した大きさ」を見ます。違いは、PCAの w がデータから学習されるのに対し、waveletは母waveletを事前に定義することです。
「異なる周波数なら反応ゼロ」ではない
無限長の理想的なsin/cosなら異周波数は直交します。しかし実際のwaveletは有限時間なので、近い周波数にも反応します。
だから、スカログラムの1点を「完全に1つのHz」と読まない。時間・周波数に幅を持つ局所的な測定です。
覚え方：W(a,b)は「時刻bの一点の値」ではありません。b周辺の窓で、scale aに対応する形がどれくらい信号に含まれるかです。
CHAPTER 4 · COMPLEX COEFFICIENT
5. 複素係数からpowerとphaseを分ける
複素Morletはcosとsinの両方を使うため、同じ周波数でも周期のどの位置にいるかを表現できます。
ψ(t) ≈ Gaussian(t) × eiωt = Gaussian(t) × [cos(ωt) + i sin(ωt)]
|W|
振幅的な強さ
局所的な成分の大きさ。スケールや正規化の影響を受けます。
|W|²
wavelet power
スカログラムの色に使うことが多い量。生理学的な重要度そのものではありません。
arg W
phase
時刻b・scale a周辺の波が周期のどこにいるか。powerが弱い領域では不安定です。
位相差は円環データ
- 0°付近：同位相
- 90°付近：1/4周期ずれ
- 180°付近：逆位相
- 179° と −179°：算術平均は0°ではない。円環上では近い。
powerは強さ、phaseは周期上の位置。位相差の平均には円統計を使う理由を説明できる。
CHAPTER 5 · CWT / DWT
6. CWTとDWTを使い分ける
どちらが優れているかではなく、欲しい出力が「細かな時間–周波数地図」か「効率的な階層分解」かで選びます。
CWT：Continuous Wavelet Transform
- scaleと時間位置を細かく走査
- 時間×周波数の連続的なマップ
- 冗長だが局所構造が見やすい
- power、phase、coherenceに向く
DWT：Discrete Wavelet Transform
- scaleと位置を離散的に間引く
- detail / approximationへ階層分解
- 圧縮、denoising、特徴抽出に向く
- CWTより局所マップは粗い
歩行研究での入口：「歩行周期のどの時期に、どの帯域の活動が変わるか」を見たいなら、まずCWTとスカログラムが直感的です。CWTを選んだ後に、統計設計を別途決めます。
CHAPTER 6 · SCALOGRAM
7. スカログラムを「地図」として読む
横軸＝時間、縦軸＝周波数、色＝power。図本体・図の構成要素・脚注を分け、赤い引き出し線なしで地図を読めるようにします。
図中ラベル：持続帯図中ラベル：ridge図中ラベル：島図中ラベル：COI
注（脚注）：高powerは「その時刻・帯域でwaveletテンプレートとの一致が強い」という意味です。神経制御、複雑性、因果を単独で示しません。脚注は図本体の構成要素ではなく、解釈上の注意です。
地図を読む順番
- まず軸と単位、周波数の上下を確認。
- 図の構成要素を「持続帯・ridge・島・COI」として読む。
- COIとpowerの弱い領域を除外候補として印を付ける。
- 元波形・イベントマーカー・測定系と照合する。
色から飛躍しない
図で見つける → 仮説を作る → 独立データまたは事前ROIで検証する、の順を守ります。
「赤いから重要」ではなく、「軸・単位・局所power・COI」を先に確認します。
図本体・図の構成要素・脚注を区別し、持続帯・ridge・島・COIを自分の言葉で説明できる。
CHAPTER 7 · WORKFLOW
8. 実践フロー：CWTは統計ではない
CWTは信号を時間–周波数表現へ変換する方法です。仮説、特徴量、統計、解釈を分離すると、解析が迷走しにくくなります。
問いを固定
何が、いつ、どの帯域で変わる？ EMG power、運動学power、phase、coherence?
データ準備
sampling、欠損、filter、artifactを確認。見たい帯域を消さない。
CWT
母waveletと周波数範囲を決め、W(f,t)を算出。COIも記録。
特徴量
power、phase、phase difference、coherence、ROI平均など。
正規化
絶対powerかrelative powerか。被験者間の振幅差をどう扱うか。
被験者単位化
複数strideは平均・頑健集約、または階層モデルへ。
統計
事前ROIなら1変量。全マップならcluster permutation等。
統合解釈
力学・臨床・artifact・再現性と合わせて生理仮説を評価。
最重要の分岐：スカログラムを見てから「一番違う場所」をROIにし、同じデータで検定しない。事前ROI、独立データ、または全マップを対象にした多重比較補正を使います。
CHAPTER 8 · PYTHON
9. Python実装：まずは1本の信号で確認
コードを動かす前に、sampling frequency、wavelet、周波数範囲、COI、正規化方針をメモします。PyWaveletsではsampling periodを渡すと周波数をHzで扱えます。
A. PyWavelets：CWT → power → phase
import numpy as np
import matplotlib.pyplot as plt
import pywt
fs = 200
dt = 1 / fs
t = np.arange(0, 10, dt)
rng = np.random.default_rng(42)
# 3 Hzは全時間、12 Hzは5秒付近だけ存在する人工信号
signal = np.sin(2 * np.pi * 3 * t)
envelope = np.exp(-0.5 * ((t - 5.0) / 0.6) ** 2)
signal += 0.8 * envelope * np.sin(2 * np.pi * 12 * t)
signal += 0.15 * rng.standard_normal(len(t))
wavelet = "cmor1.5-1.0"
target_freqs = np.geomspace(1, 30, 80)
scales = pywt.frequency2scale(wavelet, target_freqs / fs)
coef, freqs = pywt.cwt(
signal, scales, wavelet,
sampling_period=dt, method="fft", precision=12,
)
power = np.abs(coef) ** 2
phase = np.angle(coef)
plt.pcolormesh(t, freqs, power, shading="auto")
plt.yscale("log")
plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.title("Scalogram: wavelet power")
plt.colorbar(label="Power")
plt.show()最初の検証では、人工信号の「3 Hzの帯」と「5秒付近の12 Hzの島」が見えるか確認します。実データへ進む前に、周波数軸・単位・filter・COIをチェックします。
B. MNE-Python：全time–frequency mapを比較する入口
import numpy as np
from mne.stats import permutation_cluster_test, combine_adjacency
# group_A_power: (n_subjects, n_freq, n_time)
# group_B_power: (n_subjects, n_freq, n_time)
n_freq, n_time = group_A_power.shape[1:]
adjacency = combine_adjacency(n_freq, n_time)
F_obs, clusters, cluster_p, H0 = permutation_cluster_test(
[group_A_power, group_B_power],
n_permutations=5000,
threshold=None,
tail=1,
adjacency=adjacency,
out_type="mask",
seed=42,
n_jobs=-1,
)
significant_mask = np.zeros((n_freq, n_time), dtype=bool)
for cluster, p in zip(clusters, cluster_p):
if p < 0.05:
significant_mask |= clusterこのコードは統計設計の完成形ではありません。独立単位、群デザイン、閾値、置換回数、effect size、COI除外を研究計画に合わせて決めます。
CHAPTER 9 · WALKING RESEARCH
10. 1歩行周期の6チャンネルを実際に解析する
初学者が解析の流れを追えるように、3筋のEMGと3つの関節運動学を1歩行周期（約0.83秒）だけ人工的に作り、歩行phase 0–100%で各チャンネルのMorlet型CWTを計算します。図の色は計算した局所powerです。
SITUATION 1 · 今どういうデータを解析するか
1歩行周期の6チャンネル
- fs：200 Hz
- 記録：1歩行周期・約0.83秒・167点／チャンネル
- 表示軸：gait phase 0–100%
- 値は説明用の任意単位。ノイズも少し加えています。
筋活動のダミー
- 大殿筋：35 Hzをphase 55–75%にburst
- 外側広筋：45 Hzをphase 5–22%にburst
- TA：35 Hzをphase 65–85%にburst、0–15%にも弱いburst
関節角速度のダミー
- 股関節：1周期の滑らかな形＋phase 42–62%の8 Hzイベント
- 膝関節：1周期の形＋phase 25–45%の6 Hz、55–70%の12 Hzイベント
- 足関節：1周期の形＋phase 82–96%の8 Hz、30–45%の18 Hz帯
x(t) → W(f,t) → power(f,t)=Re(W)2+Im(W)2 → 表示値 log10(power+ε)
この図の位置づけ：1周期に縮めることで、各筋のburst位置と関節運動学の局所イベントを見比べやすくしています。ここでは答えを知ったうえで信号を設計し、計算結果がその特徴を回収できるかを確認します。注意：1周期だけでは1.2 Hzのような低周波を安定して推定できません。実戦では複数strideを記録し、周波数解析と周期正規化を目的別に行います。
SITUATION 2 · 結果の図
6チャンネル：生波形 → CWT power map
各パネルは上段が入力した生波形、下段が周波数×歩行phaseのpower mapです。横軸は0–100%の1周期、計算は各チャンネル内でlog-powerを表示範囲に正規化しています。
大殿筋計算中…
外側広筋計算中…
TA計算中…
股関節計算中…
膝関節計算中…
足関節計算中…
6チャンネルのCWTを計算しています…
SITUATION 3 · 図の解釈
phase 55–75%の高周波帯
35 Hz付近の明るい帯が、設計したburstの歩行周期位置に並びます。まずは「この時間帯の局所powerが高い」と記述し、筋の因果や機能を断定しません。
phase 5–22%のburst
45 Hz付近の帯がcycleの前半に反復します。大殿筋やTAと同じHzを探すのではなく、各筋のburst時期と帯域を別々に確認します。
65–85%＋弱い0–15%
TAは主burstと弱いburstが別の周期位置に現れます。raw EMG、電極配置、動作アーチファクトを併せて本物の活動かを確認します。
滑らかな形＋phase 42–62%の8 Hz
1周期の角速度の形に加えて、中央付近へ局所化した8 Hz成分を入れています。まずは「どのphase・帯域に局所powerがあるか」を読みます。
phase 25–45%の6 Hz＋55–70%の12 Hz
1周期の形に2つの局所イベントを入れています。イベントの意味は、角度・角速度・外部イベントのphaseと照合します。
phase 82–96%の8 Hz帯
終盤に局所化した帯があり、周期位置としては足部の特定局面に対応します。ただし、微分ノイズやheel-strikeの影響を除外してから解釈します。
1周期の図をどう読むか
今回の横軸：0–100%のgait phase。筋burstや関節イベントの「周期内の位置」を比べます。
実戦：実時間の秒・Hzを保った解析と、複数strideを0–100%へそろえる解析を目的別に使い分けます。
1周期の図は位置を理解する教材です。低周波の安定した推定や群比較には、複数周期と被験者単位の集約が必要です。
この図の色の注意
各パネルはチャンネル内のlog-powerを見やすく正規化しています。したがって、大殿筋の赤と膝関節の赤を「絶対powerが同じ」と比較しません。
群比較では、被験者単位の特徴量、同じ前処理、ROI、COI除外、多重比較を事前に決めます。
| 今回の表示 | 実戦で増やすもの | 向いている見方 | 注意 |
|---|---|---|---|
| 1周期・0–100% phase | 複数stride | burstやイベントの周期内位置 | 低周波の推定や群比較には不足 |
| 実時間（秒） | 連続記録・複数周期 | 物理Hz、impact、イベント時刻 | 歩行速度やstride時間の違いを別途扱う |
| 両方を保存 | cycle正規化版も作る | 物理Hzとgait phaseを分けて報告 | 変換手順・ROI・COIをMethodsに書く |
この結果から言える
「この人工1周期では、大殿筋・外側広筋・TAの高周波powerと、3関節の局所イベントが設計したphase位置に現れた」
この結果だけでは言えない
「実際の患者で同じ筋活動が起きる」「power増加が歩行改善の原因である」「6チャンネル間のpowerの大きさをそのまま比較できる」
SITUATION 4 · 筋同士の位相差を読む
大殿筋とTAは35 Hz帯でどう協調するか？
同じ試行の大殿筋とTAについて、両方に共通する35 Hz帯を選び、歩行phase 65–75%で「波の周期上の位置」がどれだけずれるかを調べます。
今回の人工データには、TA側へ約45°の位相オフセットを入れています。計算がその設計を回収できるかを確認します。
同じ周波数・十分なpower・円統計
- 同じ周波数の複素CWT係数を使う。
- 両筋のpowerが低い時刻とCOIを除外する。
- 位相差を角度として円平均し、平均角と集中度Rを報告する。
Δφ(f,t)=arg(WTA(f,t) · W大殿筋*(f,t)) R=|ΣeiΔφ|／n
実計算：35 Hzの筋間位相差
左は1歩行周期内の位相差、右はROI内の角度分布です。位相差はTA − 大殿筋として表示します。
平均位相差計算中…
集中度 R計算中…
採用サンプル計算中…
位相差を計算しています…
同じ35 Hz成分
2筋に共通する帯域だけを比較します。45 Hzの外側広筋をこの位相差へ混ぜません。
線と円で確認
phase 65–75%の線が一定の角度付近に集まり、円環上の点も一方向へ集まれば、ROI内の位相関係が安定しています。
局所的な時間関係
「35 Hz・phase 65–75%でTAが大殿筋より約○°先行／遅延」と記述します。これは同期の手がかりです。
因果ではない
位相差だけで筋が筋を駆動したとは言えません。cross-talk、power閾値、フィルタ、イベント同期を確認します。
実戦での判断：筋間位相差は解析できます。ただし「位相差がある」だけでなく、同じ周波数を比較したか、両信号に十分なpowerがあるか、位相差の分布が安定しているか、事前にROIを決めたかをセットで報告します。複数stride・被験者を比べるときは、各被験者の円統計量を作ってから群比較します。
SITUATION 5 · 論文型ケースレポート：複数strideで集約する
論文の方法を簡略化して再現
200 Hz · 12 strides / condition1周期の教材を、12 strideの研究デザインへ広げる
1周期図で「どこに帯があるか」を理解したあと、研究では通常、連続歩行から複数strideを切り出して1 strideごとの指標にします。ここでは2つの人工条件（協調良好様／協調低下様）を各12 stride作り、論文型のEMG指標を3筋で計算します。
活動が何strideで現れたか
筋ごとに仮説帯域とgait-phase ROIを決め、各strideでpowerが閾値を超えたかを判定します。出現頻度（OF）＝出現stride数／全stride数 × 100です。
OF = # {stride : peak power ≥ threshold} / Nstride × 100
活動時の周波数は何Hzか
出現したstrideのROI内でpowerを周波数方向に重み付けし、代表周波数を求めます。ここでは説明しやすいpower-weighted mean frequencyを使います。
f̄ = Σ f · power(f,t) / Σ power(f,t)
2筋の同期量を帯域で比較
大殿筋とTAの複素係数からcross-waveletを作り、30–60 Hzのcross-wavelet powerを5–80 Hz全体で割ります。論文でいうPiper帯の考え方です。
WXY = WX · WY* S30–60 = Σ|WXY|30–60 / Σ|WXY|5–80 × 100
1
入力：連続歩行をstrideへ分割
2条件 × 12 stride。各strideのphaseを0–100%にそろえます。
2
計算：CWTをstrideごとに集約
EMGの仮説帯域、power閾値、cross-wavelet帯域を先に固定します。
3
解釈：被験者単位で比較
strideは反復測定。研究では被験者内で集約してから群・条件を比較します。
図A · CWTから得た筋活動指標
出現頻度（%）と周波数内容（Hz）
図B · 2筋の同期量
normalized cross-wavelet powerの帯域構成
出現頻度（OF）計算中…
周波数内容計算中…
30–60 Hz同期量計算中…
OFが下がる＝毎strideの再現性が下がる
計算結果を読み込んでいます。
同期量は帯域内の共同活動
計算結果を読み込んでいます。
論文型の複数stride計算を実行しています…
論文との対応と差分：Di Nardo et al.の歩行sEMG研究は、TA・外側腓腹筋・大腿直筋・大腿二頭筋を複数strideでCWTし、筋活動の出現頻度と周波数内容を被験者・群で比較しました。Lodha et al.の歩行同期研究は、EMGをdetrend後に5 Hz high-passし、Morlet ω₀=6で各筋を変換、WXY=WXWY*から5–300 Hzの帯域別normalized cross-wavelet powerを計算しています。このデモは手順を追う人工データで、前処理・被験者統計は省略し、fs＝200 Hzのため5–80 Hzへ縮約しています。
方法の参照：Di Nardo et al.（歩行sEMGのCWT：出現頻度・周波数内容） ／ Lodha et al.（歩行中のMorlet cross-wavelet：30–60 Hz Piper synchrony）。論文のサンプル数・結果を再現するものではなく、解析手順を追うためのケースです。