第4章音響信号の等価性判定とピアソン相関係数¶

本章では、2つの有限長さ離散音響信号が物理的および聴覚的に実質同一であるかを数学的に判定するための基礎理論を展開する。また、計算機上で数式を評価する際に不可避となる浮動小数点演算の桁落ちを回避するためのオンラインアルゴリズムの数学的証明と、SIMD命令を用いた実装極限の最適化技術について解説する。

Note

本章における主要な略称・用語

SIMD: Single Instruction, Multiple Data
AVX: Advanced Vector Extensions
CLR: Common Language Runtime (.NET 実行環境)
FMA: Fused Multiply-Add
スケール・シフト不変性: 線形変換に対する不変性 (Scale and Shift Invariance)
DCオフセット: 直流成分のズレ

4.1 離散時間信号の類似度空間の定式化¶

音響信号の比較を行うための土台として、信号を内積が定義された計量ベクトル空間（ヒルベルト空間の有限次元の部分空間）上の点として定義する。

定義 4.1 (離散時間音響信号ベクトル)¶

サンプリング周波数 $f_s \in \mathbb{R}^+$（例えば $44100\text{ Hz}$）において、量子化された離散時間音響信号を $N$ 次元の実数ベクトル空間 $\mathbb{R}^N$ 上の点 $\mathbf{x} = [x_0, x_1, \dots, x_{N-1}]^T$ として定義する。ここで各要素 $x_i \in \mathbb{R}$ は時間ステップ $t_i = i / f_s$ における振幅を表す。実装上は $[-1.0, 1.0]$ の範囲に正規化された単精度浮動小数点数（float）として扱われる。

定義 4.2 (ピアソン積率相関係数)¶

長さ $N$ の2つの離散信号ベクトル $\mathbf{x}, \mathbf{y} \in \mathbb{R}^N$ において、その標本平均をそれぞれ $\bar{x} = \frac{1}{N}\sum_{i=0}^{N-1}x_i, \bar{y} = \frac{1}{N}\sum_{i=0}^{N-1}y_i$ とする。このとき、両ベクトル間のピアソン積率相関係数 $r_{xy} \in [-1.0, 1.0]$ は以下のように定義される。

\[r_{xy} = \frac{\sum_{i=0}^{N-1} (x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=0}^{N-1} (x_i - \bar{x})^2} \sqrt{\sum_{i=0}^{N-1} (y_i - \bar{y})^2}} \quad \cdots (4.1)\]

空間 $\mathbb{R}^N$ における標準内積 $\langle \mathbf{u}, \mathbf{v} \rangle = \sum_{i=0}^{N-1} u_i v_i$ およびそこから誘導される $L^2$ ノルム $\|\mathbf{u}\|_2 = \sqrt{\langle \mathbf{u}, \mathbf{u} \rangle}$ を導入する。平均減算した中心化ベクトル $\tilde{\mathbf{x}} = \mathbf{x} - \bar{x}\mathbf{1}$、$\tilde{\mathbf{y}} = \mathbf{y} - \bar{y}\mathbf{1}$ （ここで $\mathbf{1}$ は全要素が $1$ のベクトル）を用いると、式 (4.1) の分子は $\langle\tilde{\mathbf{x}}, \tilde{\mathbf{y}}\rangle$、分母は $\|\tilde{\mathbf{x}}\|_2 \|\tilde{\mathbf{y}}\|_2$ と書き換えられる。

したがって、相関係数は $\mathbb{R}^N$ 空間における2つの中心化ベクトルのなす角 $\theta$ のコサイン類似度に一致する。

\[r_{xy} = \frac{\langle\tilde{\mathbf{x}}, \tilde{\mathbf{y}}\rangle}{\|\tilde{\mathbf{x}}\|_2 \|\tilde{\mathbf{y}}\|_2} = \cos \theta \quad \cdots (4.2)\]

判定対象となる2つの信号の相関係数 $r_{xy}$ が、あらかじめ定めた許容閾値 $1 - \epsilon$ 以上であれば、両者は実質的に同値（線形従属に近い状態）であるとみなす。

4.1.1 ピアソン相関係数のスケール・シフト不変性と正規化パスの削減¶

前述の式 (5.1) において、ピアソン相関係数が持つ極めて強力な数学的性質としてスケール・シフト不変性が挙げられる。

任意の定数 $a > 0$ および $b \in \mathbb{R}$ に対して、信号 $\mathbf{y}$ を線形変換した $\mathbf{y}' = a\mathbf{y} + b\mathbf{1}$ を考える。このとき、変換後の平均は $\bar{y}' = a\bar{y} + b$ となり、各要素の偏差は $y'_i - \bar{y}' = (ay_i + b) - (a\bar{y} + b) = a(y_i - \bar{y})$ となる。これを式 (4.1) に代入すると、分子からは $a$ がくくり出され、分母の各項からは $\sqrt{a^2} = a$ がくくり出されるため、約分されて結果は完全に一致する。

\[r_{x, ay+b} = r_{xy} \quad (a > 0)\]

【計算量最適化への応用】 この数理的性質は、オーディオ信号処理において決定的な意味を持つ。比較する2つのWAVファイル間で、音量差やDCオフセットが存在していたとしても、波形の形状さえ一致していれば相関係数は完全に等しくなる。これにより、素朴な実装で不可避となる「比較前に全音声ファイルを一度スキャンしてピーク音量を探索し、DCオフセットを除去する」という前処理パスを丸ごとスキップすることが可能となり、ファイル読み込みから比較までの物理的なメモリアクセス量を大幅に削減している。

4.2 計算量と精度のジレンマを解消するオンライン累積法¶

4.2.1 課題：なぜ素朴な計算アルゴリズムは破綻するのか¶

式 (4.1) をプログラム上で評価する場合、計算量と計算精度の相反するジレンマに直面する。

1. 定義通りの計算（2パス・アルゴリズム）: 式 (4.1) をそのまま実装すると、まず平均 $\bar{x}, \bar{y}$ を計算するために1回目の配列走査を行い、次に偏差の積和を計算するために2回目の走査が必要となる。時間計算量こそ $O(N)$ であるが、巨大な音声データをメモリからCPUキャッシュへ2回ロードする必要があり、メモリアクセスのボトルネック（定数倍の悪化）によりパフォーマンスが著しく低下する。

2. 展開式による1パス化と桁落ち: メモリアクセスを1回（1パス）に削減するために、分散の定義式を展開した以下の等式が広く知られている。

\[\sum_{i=0}^{N-1} (x_i - \bar{x})^2 = \sum_{i=0}^{N-1} x_i^2 - \frac{1}{N}\left(\sum_{i=0}^{N-1} x_i\right)^2 \quad \cdots (4.3)\]

この式を用いれば、時間計算量 $O(N)$ かつ配列走査1回で済む。しかし、浮動小数点演算（IEEE 754）においてこの式を評価すると、桁落ちという致命的な精度崩壊を引き起こす。波形の振幅変動が極めて小さい場合、右辺の第一項と第二項は巨大な近似値となる。この減算において有効ビットの大半が相殺され、最悪の場合は丸め誤差により分散がマイナスになる。

4.2.2 解決策としての Welford の漸化式 (Why)¶

上記のジレンマを打破し、1パスによるメモリアクセス半減と桁落ちを回避する数学的安定性を両立させるため、本システムでは式を Welford のアルゴリズムに準ずる漸化式へと変形する。

この変形の本質は、配列全体をメモリに保持する（空間計算量 $O(N)$）のではなく、直前までの状態（平均・分散）と現在の入力値の微小な差分のみを用いて状態を更新していく点にある。これにより空間計算量を $O(1)$ に圧縮し、時間計算量 $O(N)$ の厳密な1パスを実現する。

4.2.3 定理 4.1 とその証明¶

Important

定理 4.1 系列の $k$ 番目までの標本平均を $\bar{x}_k$、偏差平方和を $S_{xx, k}$ とする。初期条件を $\bar{x}_0 = 0, S_{xx, 0} = 0$ としたとき、$k \ge 1$ について次式が成立する。

平均の漸化式: $$\bar{x}_k = \bar{x}_{k-1} + \frac{x_k - \bar{x}_{k-1}}{k} \quad \cdots (4.4)$$
偏差平方和の漸化式: $$S_{xx, k} = S_{xx, k-1} + (x_k - \bar{x}_{k-1})(x_k - \bar{x}_k) \quad \cdots (4.5)$$

Note

証明

平均の漸化式について：定義より $\bar{x}_k = \frac{1}{k}\sum_{i=1}^k x_i$ である。これを変形すると、

\[ \begin{aligned} \bar{x}_k &= \frac{1}{k} \left( x_k + \sum_{i=1}^{k-1} x*i \right) \\ &= \frac{1}{k} \left( x_k + (k-1)\bar{x}*{k-1} \right) \\ &= \frac{x*k + k\bar{x}*{k-1} - \bar{x}_{k-1}}{k} \\ &= \bar{x}_{k-1} + \frac{x*k - \bar{x}*{k-1}}{k} \end{aligned} \]

となり、(4.4) が示された。
偏差平方和の漸化式について： $S_{xx, k} = \sum_{i=1}^k (x_i - \bar{x}_k)^2$ である。和を $k-1$ 番目までと $k$ 番目に分割し、$\bar{x}_k$ を $\bar{x}_{k-1}$ を用いて表す過程で (4.4) より導かれる関係式 $(x_i - \bar{x}_k) = (x_i - \bar{x}_{k-1}) - \frac{x_k - \bar{x}_{k-1}}{k}$ を用い、展開と整理を繰り返すことで、以下の簡潔な差分更新式が得られる。 $$S_{xx, k} - S_{xx, k-1} = (x_k - \bar{x}_k)(x_k - \bar{x}_{k-1})$$ これにより (4.5) が示された。共分散 $S_{xy, k}$ についても同様に $S_{xy, k} = S_{xy, k-1} + (x_k - \bar{x}_{k-1})(y_k - \bar{y}_k)$ が導かれる。$\blacksquare$

4.2.4 実装への射影：オーダー削減と疑似コード¶

C#における PearsonAccumulator は、この定理を定数空間 $O(1)$ で実行するアキュムレータとして実装されている。

// 【疑似コード】PearsonAccumulator による O(N)時間・O(1)空間アルゴリズム
double meanX = 0, meanY = 0;
double varX = 0, varY = 0, covXY = 0;

// キャッシュラインに乗った配列を1度だけ走査する（1パスによるメモリ帯域幅の半減）
for (int i = 1; i <= N; i++) {
    double x = signalX[i];
    double y = signalY[i];

    double dx = x - meanX; // 式(4.4)の分子
    double dy = y - meanY;

    meanX += dx / i;       // 平均の更新 O(1)
    meanY += dy / i;

    varX += dx * (x - meanX); // 式(4.5)の適用 O(1)
    varY += dy * (y - meanY);
    covXY += dx * (y - meanY);
}
// 最終的な相関係数の導出
double r = covXY / Math.Sqrt(varX * varY);

この実装により、一時配列のメモリアロケーションが完全に排除され（空間 $O(N) \to O(1)$）、かつ $(x_k - \bar{x}_{k-1})$ という「真値に近い微小な差分」を累積するため、式(4.3)のような巨大数の減算が発生せず、二重の課題（速度と精度）が解決されている。

4.3 時間シフトを伴う相互相関計算とハードウェア最適化¶

実際の音響ファイル同士を比較する際、サンプリングの開始位置が微小に異なる場合がある。このとき、単なる要素ごとの比較（ラグ $\tau = 0$）では相関係数が著しく低下する。

定義 4.3 (ラグ $\tau$ を伴う相互相関関数)¶

ベクトル $\mathbf{y}$ を時間軸方向に $\tau \in \mathbb{Z}$ サンプルだけシフトさせた信号ベクトルを考える。比較可能なオーバーラップ区間長を $M(\tau) = N - |\tau|$ とする。この有限区間におけるピアソン相関係数を $r_{xy}(\tau)$ と定義する。等価性判定問題は、許容される最大のズレ幅を $\tau_{\text{max}}$ としたとき、探索範囲 $\tau \in [-\tau_{\text{max}}, \tau_{\text{max}}]$ における最大相関を見つける局所最適化問題として定式化される。

\[r_{\text{max}} = \max_{-\tau_{\text{max}} \le \tau \le \tau_{\text{max}}} r_{xy}(\tau) \quad \cdots (4.6)\]

4.3.1 【実装への射影: `FastWaveCompare`】SIMDによる計算量定数倍の極限削減¶

ラグ $\tau$ の全探索は、素朴な実装では時間計算量 $O(N \cdot \tau_{\text{max}})$ を要求する。これを実用的な時間内に収めるため、FastWaveCompare.cs ではアルゴリズム的オーダー削減に加え、ハードウェアレベルでの定数倍の極限削減を行っている。

ポインタによる境界チェック排除 C# の unsafe と fixed により配列への直接ポインタを取得する。これによりCLRの配列境界チェックループごとのオーバーヘッド $O(1)$ が消失する。
SIMD超並列化による命令数削減 スカラー演算による $\sum$ の評価を棄て、Vector256<float> を用いる。これにより、8次元のベクトル積和演算を1命令で実行する。ループの回転数は $N$ から $N/8$ へと削減され、計算スループットは理論上最大8倍に引き上げられる。

// 【疑似コード】AVX2 / AVX-512 による内積計算のSIMD最適化
Vector256<float> acc = Vector256<float>.Zero;
for (int i = 0; i < N - 7; i += 8) {
    // 8サンプルを同時にレジスタへロードし、1クロックで積和(FMA)演算
    Vector256<float> vX = Avx.LoadVector256(pX + i);
    Vector256<float> vY = Avx.LoadVector256(pY + i + tau);
    acc = Avx.Add(acc, Avx.Multiply(vX, vY));
}

ループアンロール さらに、GenerateSimdBatchUnrollAttribute を通じてループ展開を強制する。これにより、複数のSIMDレジスタに対する独立した命令がパイプラインに投入され、CPUのアウト・オブ・オーダー実行が誘発される。依存関係によるストール（待機状態）が排除され、レイテンシという物理的な壁を越えた最適化が達成される。

4.4 等価性判定と統計的確実性限界¶

算出された最大相関係数 $r_{\text{max}}$ から最終的な同一性判定を下す際、単なる固定の閾値（例えば $r > 0.99$）を用いることは統計学的に不十分である。

WaveValidation.Pearson.cs では、オーバーラップした有効サンプル数 $M(\tau)$ が少ない場合、無相関のホワイトノイズ同士であっても偶然に高い相関が生じる確率（偽陽性率）が上昇するという数理統計学的性質を考慮している。そのため、判定関数 IsPearsonCorrelationValid では、経験的かつ数学的限界に基づく厳格な許容誤差（$\epsilon \approx 10^{-4}$ オーダー）の適用に加え、十分な $M(\tau)$ の長さを必須条件として組み込むことで、確実な等価性証明を担保している。

4.5 参考文献¶

Welford, B. P. (1962). "Note on a method for calculating corrected sums of squares and products". Technometrics, 4(3), 419-420.

第4章 音響信号の等価性判定とピアソン相関係数¶