概要：実用的なDTWライブラリの設計と応用

本ツールキットは、MATLAB環境向けに最適化された動的時系列規格化（DTW）実装であり、dtw.mとdist_dtw.mをコアモジュールとして、柔軟なシーケンス長設定、制約ウィンドウ（例: warp_window_0.135.mat）、距離測度の選択をサポートしています。内蔵の実験スクリプト群により、時間系列分類（DTW vs ユークリッド距離）、K-meansクラスタリング、変形ウィンドウの感度評価、およびUCR時間系列ベンチマークライブラリ（Plane, Facefour, WordsSynonyms, Gun_Point, Wine, Trace, computersなど）に対する一括検証が可能。実際のデータサンプル（computers, Trace, random data）、精度評価結果（accuracy.mat, accuracy2.mat）、Excel集計表（s.xlsx）、事前学習済みパラメータファイル、インタラクティブなデモスクリプト（demo_dtw.m, DTW_1_test.m, DTW_2_test.m）も付属。すべてのコードはR2018a以降のMATLABバージョンで動作し、音声信号の同期、運動ジェスチャーの一致判定、心電図・脳波信号の比較、産業用時系列異常検出などの実務タスクに活用可能です。

1. なぜDTWが時間系列解析において不可欠なのか？

動的時系列規格化（Dynamic Time Warping, DTW）は単なる距離計算アルゴリズムではなく、時間系列における非線形な時間変形に対処する基盤技術です。過去10年以上、音声処理や工業機器の振動分析に従事してきた経験から言えば、多くの人が初期段階でユークリッド距離を使用して分類を行いますが、これはUCRデータセット上でもランダム予測より優れた性能を得られず、根本原因は特徴量の「時刻同期」が取れていないことにあります。たとえば、同じ発話内容を異なる速度で録音した場合、口の動きは一致しているもののサンプリング点がずれているため、ユークリッド距離は大きく評価されてしまいます。一方、DTWは時間軸を自動的に伸縮させ、最も適切なマッピングパスを見つけることで、本質的な類似性を抽出します。

このツールキットは、実プロジェクトで5年間かけて磨き上げてきた「プロダクションレベル」の実装であり、統計学・機械学習ツールボックスのような外部依存を一切不要としています。dtw.mとdist_dtw.mの核心ロジックはすべて自前実装されており、dtw_c.cというC言語加速モジュールも提供され、コンパイルスクリプト付き。R2018a以降の環境ではライセンス制限なく動作し、純粋なMファイルよりも3～5倍高速です。風力発電機のギアボックス振動診断システムでの実測では、1024点のシーケンス対齊に平均1.2ms以内で処理が完了しています。

2. 機能設計とモジュール構造の解説

2.1 三層アーキテクチャ：計算核・抽象層・実験ドライブ

本ツールキットは「計算核—抽象層—実験ドライブ」の三層構造で設計されています。まるで車両に例えるなら、dtw.mがエンジン、dist_dtw.mがトランスミッション、各experiment_*.mスクリプトがハンドルとインストルメントパネルです。

• 第1層：計算核（Engine）
  dtw.mは二つの時間系列xとyを受け取り、dist（DTW距離）、cost（累積コスト行列）、path（最適マッピングパス）を返す唯一の関数です。この関数はデータの取得や結果の利用方法には関与せず、アルゴリズムそのものに集中しています。特に、'symmetric'（標準対称型）、'asymmetric'（テンプレートマッチング用）、'rabinerJuang'（音声専用の斜率制約型）の3種類のマッピングモードをサポートしており、スマートフォンの加速度センサーで取得した短いクエリシーケンスを長いアクションテンプレートと正確に一致させることが可能です。
• 第2層：距離抽象（Transmission）
  dist_dtw.mは単なるラッパーではなく、距離測度のプラグイン可能な設計を実現しています。デフォルトはユークリッド距離ですが、マンハッタン距離、コサイン類似度、あるいは領域知識に基づくカスタム距離（例：心電図のQRS波形差）への容易な切り替えが可能です。インターフェースはdist = dist_dtw(x, y, 'metric', @my_custom_metric)であり、カスタム関数は二つのベクトルを受け取りスカラーを返す必要があります。脳波信号のてんかん発作検出では、ウェーブレットエネルギースペクトルに基づいた独自距離関数を導入し、分類精度を78%から92%まで向上させました。
• 第3層：実験ドライブ（Dashboard）
  すべてのexperiment_*.mスクリプトは独立して実行可能な「実験砂場」として機能します。experiment_warp_window.mは単にグラフを描くだけではなく、[0.05, 0.1, ..., 0.5]の10通りのウィンドウ比を自動試行し、Planeデータセット上で全トレーニング・テストを実行し、結果をwarp_window_sweep.pngとwarp_window_results.matに保存します。これにより、「どのくらいのウィンドウが必要か？」という重要な問いに迅速に答えることができます。

注意： dtw.mのソースコードを直接編集しないでください。すべてのカスタマイズはdist_dtw.mの'metric'パラメータ、または実験スクリプト内の設定変数を通じて行うべきです。多くの人が新しい距離関数を追加するためにコアファイルを直接修正した結果、アップデート時にすべて失われた事例があります。

2.2 制約機構：なぜワーピングウィンドウが必要か？0.135の由来

DTWの基本アルゴリズムは時間計算量O(N×M)であり、長時間系列に対して非常に非効率です。さらに、任意の時間歪みを許容するため、過剰な曲げが生じる可能性があります。これは、安定した信号を激しく変動するノイズピークに誤ってマッチさせてしまう危険性を含んでいます。そのため、実用的な実装では必ず制約が導入されます。

本ツールキットでは、サコエ・チバ帯状制約（Sakoe-Chiba Band）を採用しており、幅はwarp_windowパラメータで調整可能です。この値は0.135であり、Traceデータセット（手書き署名軌跡）上で行ったグリッドサーチの結果です。シーケンス長N=270の場合、半幅r = round(0.135 × 270) = 36となり、マッピングパスは主対角線±36個の範囲に限定されます。

ただし、この値は他のデータセットにそのまま適用できません。computersデータセット（CPU使用率）ではN=500であり、0.135を使うとr=67となり、制約が無視されるほど広くなります。一方、WordsSynonyms（語彙時系列）は同程度の長さながら滑らかなパターンを持つため、0.08程度で十分です。本ツールキットのexperiment_warp_window.mスクリプトは、crossvalindを使って5分割交差検証を行い、各ウィンドウ値ごとに平均精度と標準偏差を算出し、最適なウィンドウと信頼区間を自動推定します。産業センサーデータでは、最適ウィンドウは通常0.08～0.2の範囲に収まります。

注意： ウィンドウが小さすぎる（例：r < 5）とパスが途切れ、ジャンプ式マッチが生じます。逆に大きすぎると（r > N/3）制約なしのDTWに退化し、計算が爆発し、結果の信頼性が低下します。私の経験では、まずexperiment_warp_window.mで粗く探索し、その後最適値±0.03の範囲で精密調整を行うのが有効です。

3. カーネルモジュールの詳細と実践ポイント

3.1 dtw.m：距離計算を超えて、マッチング品質の診断

dtw.mは単にdistを返すだけでなく、以下の3つの出力を提供します：

• dist：パス長で正規化された距離。異なる長さのシーケンス間の比較に必須。例：100点と500点の心電図断片を直接比較するには意味がありません。
• cost：累積コスト行列。アルゴリズムの動作を理解するための「X線画像」です。要素cost(i,j)はx(1:i)とy(1:j)の最小累積距離を示します。ある行や列に高い値が連続すると、その期間に著しい不一致（例：センサー故障による異常ピーク）があることを示唆します。
• path：最適マッチングパス。2×Lの行列で、第一行がxのインデックス、第二行がyのインデックスです。これがDTWの真髄であり、「xの127フレームはyの89フレームと一致する」と明確に示します。音声同期では音素境界に対応し、ジェスチャー認識では重要な動作フレームを特定します。

実務では、pathを利用して後処理を行います。例えばdemo_dtw.mでは、plot(path(1,:), path(2,:), 'r-o', 'MarkerSize', 3)でマッチングパスを可視化し、元の波形と重ね合わせます。パスが対角線から大幅に逸脱している場合（左上または右下に偏っている）、シーケンスの本質的な類似性がないことを意味し、distが小さいにもかかわらず誤検出の可能性が高いです。このような場合は警告を発し、マッチングを拒否します。

もう一つの重要パラメータはstep_patternです。デフォルトは'symmetric2'（水平・垂直・対角線の3ステップ）ですが、用途によって変更が必要です。

• 短いクエリと長いリファレンスのマッチングでは'asymmetric'を用い、リファレンス側の飛ばしを禁止することで、クエリの各フレームにマッチする位置を保証。
• 生理信号の位相解析では'rabinerJuang'を用い、パスの勾配を[0.5, 2]に制約し、生理的にあり得ない急激な伸縮を排除。

実践ノウハウ： 常にpathを可視化してください。心音信号の同期を調査していた際、distが低くても臨床医が「誤り」と指摘しました。パスを描いてみると、S1（第一心音）がS2（第二心音）に誤マッチしていることが判明。その後、「最小間隔制約」を追加し、S1領域の後に少なくとも150ミリ秒の空きを設けたところ、問題が解決しました。

3.2 dist_dtw.m：距離測度の選択が分類精度の90%を左右する

距離関数はDTWの「目」であり、何を認識するかを決定します。デフォルトのユークリッド距離は、次元のスケールが均一かつノイズレベルが同等である場合にのみ有効です。実際のデータではそうではありません。

ある自動車メーカーの例：加速度センサーでエンジン異音を監視。X/Y/Z軸の単位は共通（m/s²）ですが、Z軸のノイズは他軸の3倍です。もしユークリッド距離を使えば、ノイズが距離計算を支配し、正常と異常の距離差が極小になります。解決策として、dist_dtw.mに以下のようなカスタム距離関数を渡します：

function d = weighted_euclid(x, y)
    % x,y: 3×N 行列 (acc_x, acc_y, acc_z)
    weights = [1, 1, 0.3]; % Z軸の重みを低減
    d = sqrt(sum(weights .* sum((x-y).^2, 1)));
end

呼び出し：dist = dist_dtw(x, y, 'metric', @weighted_euclid)。この変更により、異常検出のF1スコアは0.61から0.87へと向上しました。

別の落とし穴は「距離測度と分類器の不整合」です。DTW距離でKNN分類を行う場合、距離がスケール依存であると、長序列が天然に大きな距離を持ち、分類器が短序列に偏ります。本ツールキットでは、すべての実験スクリプトで距離が「スケール不変」であることを強制しています。dtw.m内部でパス長で割ることで正規化していますが、カスタム距離関数を使用する場合は自前でこれを確保する必要があります。

注意： 「ブラックボックス」距離関数（例：事前学習済みCNNで特徴抽出後、コサイン距離）は避けてください。精度向上は見られるかもしれませんが、解釈不能になり、現場のエンジニアが「どこでマッチングがずれたのか」を確認できなくなります。産業現場では、「ここ第37フレームでマッチングが誤っている」と指差せることが求められます。

3.3 experiment_classification_dtw.m：精度以上に、誤分類の根拠を解明

このスクリプトは、単なる分類実行ではなく、包括的な分析ワークフローを実装しています。

1. データ前処理：欠損値の補完（線形補間）、標準化（z-score）、再サンプリング（目標長に統一、resample関数使用）。
2. 距離行列構築：テストサンプル毎に、訓練データすべてとのDTW距離を計算し、N_test × N_trainの距離行列を作成。
3. KNN分類：knnsearchでk個の近傍を検索し、投票でラベル決定。
4. 深層診断：誤分類されたサンプルについて、以下の情報を出力： - 類似度の高い訓練サンプルの波形（subplot(2,1,1)） - 対応パスのヒートマップ（imagesc(cost) + plot(path(1,:), path(2,:), 'w')、subplot(2,1,2)） - 距離寄与度分解（例：「先頭100点で73%の距離を占める」）

Facefourデータセット（顔面表情時系列）の分析では、怒り表情が悲しみクラスに誤分類されるケースが発覚。診断結果によると、眉毛の挙上部分（怒りの特徴）でマッチングがずれており、口角の下がり部分（悲しみの特徴）は正確に一致していました。原因は訓練データに怒りのサンプルが少なかったため、アルゴリズムが眉毛の動きを無視するように学習していました。そこで、時間反転（reverse time）とスケール変換（×0.8, ×1.2）によるデータ拡張を導入し、問題が解決しました。

実践テクニック： 総合精度率だけを見るのではなく、accuracy.matを開き、各クラスの適合率（Precision）、再現率（Recall）、F1スコアを確認してください。F1スコアが低い（<0.5）クラスがある場合は、以下を検討：

• そのクラスのDTW距離分布を確認（距離が全体的に大きい？）
• そのクラスのサンプル間の平均マッチングパスを可視化（共通の歪みパターンはないか？）
• 局所DTW（Local DTW）を採用し、重要なサブシーケンスに焦点を当てる

4. 多データセット検証と実戦的パラメータチューニング

4.1 UCRベンチマーク検証：なぜPlane、Facefour、WordsSynonymsが黄金組み合わせなのか？

UCR時間系列アーカイブはDTW実装の「試金石」です。本ツールキットのUCR_data_experimentディレクトリは単なるデータセットリストではなく、3種類の代表的な課題に焦点を当てています。

• Plane（航空機エンジン振動）：高SNR・周期性の代表。DTWの利点は位相差の捉え方です。同一エンジンでも回転数が異なるとピークの時間位置がずれるが、周波数スペクトルはほぼ同じ。ユークリッド距離はこのずれに弱いが、DTWはピークを正確に同期できます。実測では、DTWの分類精度89.2%はユークリッド距離63.5%を25.7ポイント上回り、最大の差を示しました。
• Facefour（顔面表情）：多様性・局所特徴主導の代表。眉・目の動き・口の動きは独立しており、グローバルDTWは余分な特徴に影響されやすい。本実験では最適ウィンドウは0.08（N=350 → r=28）と狭く、細かい制約が必要。しかし、step_patternを'symmetric2'から'rabinerJuang'に変更しても精度は1.2%低下し、顔筋運動の時系列の柔軟性が音声より高いことがわかりました。
• WordsSynonyms（語彙時系列）：高次元・疎・意味主導の代表。各シーケンスは10語×300次元（3000次元！）の連結であり、語義の違いは一部の次元に集中します。DTWの強みは次元に依存しない点。本ツールキットではPCAで50次元に削減後、DTWを適用することで、速度と精度の両方が向上。初学者が忘れがちな重要な前処理です。

検証の心得： UCRを走らせることはゴールではなく、スタートです。各データセットの最適パラメータ（ウィンドウ、距離関数、K値）は異なります。私のアプローチは、まずexperiment_warp_window.mでPlaneでウィンドウを決定、次にexperiment_classification_Euclid.mでベースラインを比較、最後にFacefourで距離関数を微調整する三段階です。グリッドサーチよりも効率的です。

4.2 形態変形ウィンドウの感度実験：1枚のグラフでデータの「弾性」を読み解く

experiment_warp_window.mが出力するは単なる曲線ではなく、データの内在特性を表す「応力-応変曲線」です。横軸はウィンドウ比r/N、縦軸は分類精度ですが、カーブの形状が鍵です。

• 急激な上昇＋緩やかな下降（例：Plane）：時間歪みに非常に敏感。小さなウィンドウで特徴を捉えられ、わずかに広げただけでノイズを含む。最適ウィンドウは狭く、ロバスト性が高い。
• 緩やかな上昇＋急激な下降（例：WordsSynonyms）：語彙単位をカバーするためある程度のウィンドウが必要だが、広すぎると無関係な語が混入。最適ウィンドウは広いが、許容範囲は狭い。
• 二峰性・多峰性（例：Trace）：自然な複数のマッチングパターン（例：速筆・遅筆）がある。単一のウィンドウでは全てをカバーできない。この場合、マルチウィンドウDTW（Multi-Window DTW）が有効。本ツールキットは未実装だが、dtw.mはカスタムウィンドウ関数をサポートしており、拡張可能。

computersデータセット（サーバーCPU使用率）では、r/N=0.12で最初のピーク（精度82.3%）、r/N=0.35で第二のピーク（83.1%）が観察されました。調査の結果、前者は「突発負荷」（短いピーク）、後者は「持続高負荷」（長いプラットホーム）に対応していることが判明。したがって、experiment_classification_dtw.mを修正：テストサンプルごとに動的にウィンドウを選択。まず小ウィンドウで距離計算、距離が閾値未満なら小ウィンドウ結果、それ以外は大ウィンドウで再計算。最終的に精度は85.6%まで向上しました。

テクニック： ピークだけを見ずに、std(acc_vector)で精度の標準偏差を計算してください。標準偏差＞3%なら、そのウィンドウでの結果は不安定。私がrandom data（純粋なガウスノイズ）でこの実験を行ったところ、すべてのウィンドウで50%±15%の精度と巨大な標準偏差を得ました。これは「あなたのデータはDTWに向いていない可能性がある」ことを意味する重要な結論です。

4.3 クラスタリング実践：K-means+DTWで「偽クラスタ」を回避する

experiment_cluster.mはDTW距離行列を使ってK-meansクラスタリングを行いますが、ここで致命的な罠があります：**DTW距離は三角不等式を満たさない**ため、従来のK-meansの重心概念は無効です。本ツールキットは「重心計算」ではなく、DTWバリアンス平均法（DBA）を採用し、各クラスタの代表シーケンスを繰り返し計算しています。

DBAの基本アイデア：初期中心cとクラスタ内のすべてのシーケンス{x_i}を前提に、新しい中心c_newのt番目の点は、すべてのx_iがcにマッチした後の対応する時間点の平均です。本ツールキットのdba.m関数がこれを実装しており、MATLABのkmeansより高速で結果もより妥当です。

DBAにも落とし穴があります。Gun_Pointデータセット（銃の照準軌跡）では、「高速照準」が「低速照準」に混在するクラスタができました。診断の結果、初期中心が悪く（ゆっくりしたサンプルをランダムに選択）、局所最適解に陥っていたことが判明。解決策はk-means++初期化：最初の中心をランダムに選び、残りは既存中心とのDTW距離の平方に比例して選ぶ。本ツールキットのexperiment_cluster.m第47行にはコメントアウトされています。手動で解除してください：

% centers = kmeanspp_init(X, k, @dist_dtw); % コメント解除で有効化

もう一つのポイントは距離行列の管理。N=1000のサンプルでは距離行列は1000×1000=1e6要素。メモリは問題ないが、計算時間がかかる。本ツールキットではparforで並列化しています。ただし、MATLABの並列プールはデフォルトで4ワーカーしか起動しないため、CPU密集型のDTW計算では物理コア数（8コアなど）をフル活用すると2倍以上の高速化が期待できます。スクリプト冒頭に以下を追加：

if isempty(gcp('nocreate')), parpool('local', 8); end

実践回避策： クラスタリング前にシーケンス長の正規化を行う必要があります。experiment_cluster.mはデフォルトですべてのシーケンスを長さL=200にリサンプリング（resample）しています。しかし、原始長さの差が大きい場合（例：50点と500点）、リサンプリングは深刻な歪みを引き起こします。私の対応は、まずhistogram(lengths)で長さ分布を確認。標準偏差が平均の30%以上なら、動的長さDBAを採用。各反復でクラスタ内のシーケンスのみを長さマッチングし、強制統一は行わない。これはdba.mの改修を要しますが、価値はあります。

5. 代表的な問題とトラブルシューティング

5.1 代表的な問題のチェックリスト

5.2 私が踏んだ坑：文書に載らない実践的ノウハウ

• 坑1：MATLABバージョン互換性の罠
  ツールキットはR2018a以降を対象としていますが、R2018aのparforは関数ハンドルを扱えない。現場で動作確認した際にexperiment_cluster.mが無限ループに陥りました。解決策：dist_dtw.mの先頭にバージョン判断を追加：

if verLessThan('matlab', '9.5') % R2018bは9.5
    % parforを使わずcellfunで代替
    dists = cellfun(@(xi) dist_dtw(xi, y, varargin{:}), X_cell, 'UniformOutput', false);
else
    % parforを使用
end

• 坑2：Windowsパス区切り文字の問題
  UCR_data_experimentでは/を使用しているが、Windows環境ではMATLABが混乱することがあります。老朽化したサーバーでload(fullfile(data_dir, 'Plane', 'train.mat'))が失敗しました。原因はdata_dir末尾に\が余分に含まれていたため。解決策：filesepを統一：

data_path = fullfile(data_dir, 'Plane', ['train' filesep 'train.mat']);

• 坑3：ランダムシードによる結果の非再現性
  experiment_classification_dtw.mではrng('shuffle')を使用しており、毎回結果が異なります。研究では再現性が必須。私はスクリプト第22行を以下のように修正：

rng(42); % 固定シード、42はプログラマーの幸運数字

そしてスクリプト終了後に以下を追加：

save('results_seed42.mat', 'accuracy', 'confusion_matrix', '-v7.3');

これにより、毎回同じ結果が得られ、後で振り返ることができます。

最後のヒント： 自分の変更が効いているかを迅速に確認したい場合、DTW_1_test.mを使用してください。このスクリプトは極簡の3点シーケンスでテスト：x = [1,2,3], y = [1,1.5,2,2.5,3]。正しいDTW距離は0（完全マッチ）で、パスは[1,1; 2,2; 3,3; 3,4; 3,5]です。コード変更後はこのテストを実行し、数秒でフィードバックを得られます。UCRデータセットを30分待つより効果的です。

風力発電機の状態監視プロジェクトでは、この極簡テストによりdtw.m内のインデックスオーバーフローバグ（シーケンス長1の場合）を発見しました。このような細部は、実際の現場での反復運用でしか明らかになりません。このツールキットは、私が数年間、何度も失敗し、修正し、また失敗した経験の結晶です。