GCC-PHATとは？PHAT重みの意味とTDOA推定の手順・Python実装

GCC-PHAT（Generalized Cross-Correlation with Phase Transform）は、2本のマイクに届いた音の時間差（TDOA）を求めるアルゴリズムです。相互スペクトルの振幅を周波数ごとに1へそろえ、位相だけを残してから時間領域へ戻すため、相関のピークが鋭くなります。この記事では、PHATという重みの意味、計算の手順、そのまま動くNumPy実装と実行結果、SRP-PHATとの違いを順に説明します。反射音と雑音を加えたときにPHATがどう効き、どこで効かなくなるかも、実行した数値で示します。

まとめ：GCC-PHATの要点

• GCC-PHATは、一般化相互相関（GCC）の重み関数に「相互スペクトルの振幅の逆数」を使う方式です。GCCの枠組みはKnappとCarterが1976年のIEEE Transactions on Acoustics, Speech, and Signal Processing（24巻4号）で示しました。
• PHATは位相変換（Phase Transform）の略です。振幅を捨てて位相だけで相関を取るので、低域にエネルギーが偏った音声でもピークが鈍りにくくなります。
• 手順は「フレーム切り出し→FFT→相互スペクトル→振幅で割る→逆FFT→ピーク探索」の6段で、NumPyなら十数行で書けます。
• 後述の検証では、反射音ありの条件でGCC-PHATが20回中20回、重みなしの相互相関が20回中14回、真の遅れを当てました。一方、SNR 0dB以下では重みなしと同等か下回りました。
• 多数のマイクで空間を走査する場合は、GCC-PHATを全ペアで足し合わせるSRP-PHATへ進むのが定石です。

GCC-PHATの定義とPHATの意味

一般化相互相関（GCC）と重み関数

2つのマイク信号 x1、x2 の相互相関は、周波数領域では相互スペクトル G12(f) = X1(f)・X2*(f) を逆フーリエ変換したものです。KnappとCarterの論文「The generalized correlation method for estimation of time delay」（1976年、320〜327ページ）は、この相互スペクトルに重み関数 Ψ(f) を掛けてから逆変換する形を「一般化相互相関」として整理しました。重みを何にするかで、ROTH、SCOT、PHAT、ML（最尤）などの方式に分かれます。

式で書くと次のとおりです。R12(τ) が最大になる τ を、2本のマイク間の到来時間差として採用します。

R12(τ) = ∫ Ψ(f) · X1(f) · X2*(f) · e^(j2πfτ) df

PHAT（Phase Transform）重みと白色化の効果

PHATの重みは Ψ(f) = 1 / |G12(f)| です。相互スペクトルを自身の大きさで割るので、どの周波数でも大きさが1になり、残るのは位相 e^(-j2πfτ0) だけになります。理想的な単一経路なら、これを逆変換した結果は τ0 の位置に立つ鋭いインパルスになります。

重みなしの相互相関は、エネルギーの大きい周波数帯に結果が引っ張られます。音声は数百Hz以下にエネルギーが集まるため、相関のピークがなだらかに広がり、近くに反射音のピークがあると区別しにくくなります。PHATはスペクトルを平らにする白色化の一種で、電子情報通信学会の知識ベースでは同じ処理を「白色化相互相関（CSP：Cross-power Spectrum Phase Analysis）法」と呼んでいます。音響信号処理で使うPHATは、英語のスラングとは異なり、Phase Transformの略です。

画像処理の正規化相互相関（NCC）は、ウィンドウ全体のエネルギーで1回だけ割って値を-1〜1に収めます。周波数ごとに割るPHATとは、正規化する単位が違います。NCCの考え方はテンプレートマッチングとは？matchTemplateの仕組みと類似度指標の選び方を解説で扱っています。

TDOAから音の到来方向を求める計算

マイク間隔を d、音速を c、到来角（マイクを結ぶ線の垂線から測った角度）を θ とすると、遠方の音源では τ = d・sinθ / c が成り立ちます。GCC-PHATで τ が求まれば、θ = arcsin(c・τ / d) で角度に直せます。

探索すべき遅れの範囲は ±d/c に限られます。d=10cm、c=343m/sなら最大で約0.29msです。サンプリング周波数16kHzでは1サンプルが約21.4mmの経路差に当たり、±4.66サンプルしか範囲がありません。正面付近では1サンプルの違いが約12.4°の角度差になるため、サンプル単位のピーク位置だけでは粗すぎます。後述の実装で補間を入れているのはこのためです。

マイクが2本だけでは、前方から来た音と後方から来た音が同じ τ になります（前後の曖昧性）。360°の方位を一意に決めるには、同一直線上にない3本以上のマイクが必要です。マイクがM本あればペアは M(M-1)/2 通りあり、4本なら6ペアです。

GCC-PHATの計算手順（フレーミングからピーク探索まで）

1. フレーム切り出し：連続信号を数十ms程度の窓で切り出します。音源が動く場合は、窓を短くするほど追従しやすくなりますが、1フレームの情報量は減ります。
2. FFT：各マイクのフレームを複素スペクトル X1(f)、X2(f) に変換します。
3. 相互スペクトル：X1(f)・X2*(f) を周波数ごとに計算します。この位相がマイク間の遅れを表します。
4. PHAT重み：相互スペクトルを自身の絶対値で割ります。無音の帯域で0除算しないよう、分母に小さな定数を足します。
5. 逆FFT：時間領域の相関関数 R12(τ) に戻します。周波数領域で零詰めしてから逆変換すると、サンプル間の値を補間できます。
6. ピーク探索：±d/c の範囲で |R12(τ)| が最大になる τ を探し、角度へ変換します。

動く音源を追う場合、フレームごとの推定値はばらつきます。実用システムでは推定値をそのまま使わず、ヒストグラムやカルマンフィルタとは？仕組みと使い方をわかりやすく解説【線形の基礎から応用まで】のような追跡処理で平滑化します。

PythonによるGCC-PHAT実装と実行結果

NumPyだけで書くGCC-PHAT関数と検証コード

次のコードは、GCC-PHAT関数と、比較用の重みなし相互相関、後述の検証を1ファイルにまとめたものです。依存はNumPyとは？読み方・インストール・基本的な使い方をPythonで解説だけです。Python 3.9.6・NumPy 2.0.2で実行しました。

import numpy as np

def gcc_phat(x, y, fs, max_tau=None, interp=16):
    """y を基準にした x の遅れ（秒）を返す。正なら x が後から届いている。"""
    n = len(x) + len(y)
    X = np.fft.rfft(x, n=n)
    Y = np.fft.rfft(y, n=n)
    R = X * np.conj(Y)                   # 相互スペクトル
    R /= np.abs(R) + 1e-12               # PHAT重み：振幅を1にそろえる
    cc = np.fft.irfft(R, n=interp * n)   # 零詰めで補間しながら時間領域へ
    max_shift = int(interp * n / 2)
    if max_tau is not None:
        max_shift = min(int(interp * fs * max_tau), max_shift)
    cc = (np.concatenate((cc[-max_shift:], cc[:max_shift + 1]))
          if max_shift > 0 else cc[:1])
    if not np.any(cc):
        return np.nan, cc
    shift = np.argmax(np.abs(cc)) - max_shift
    return shift / float(interp * fs), cc

def plain_cc(x, y, max_lag):
    """重みなしの相互相関（探索範囲を±max_lagに制限）"""
    c = np.correlate(x, y, mode="full")
    mid = len(y) - 1
    w = c[mid - max_lag:mid + max_lag + 1]
    return np.argmax(np.abs(w)) - max_lag

fs, N, LAG = 16000, 4096, 40
TRUE = 3  # 真の遅れ（サンプル）

def source(rng):
    # 低域に偏った信号（64サンプル移動平均した白色雑音）
    return np.convolve(rng.standard_normal(N), np.ones(64) / 64, mode="same")

# 実験1：反射音（20サンプル後・振幅0.8）がある場合
p = c = 0
for s in range(20):
    src = source(np.random.default_rng(s))
    m1 = src
    m2 = np.roll(src, TRUE) + 0.8 * np.roll(src, TRUE + 20)
    p += round(gcc_phat(m2, m1, fs, max_tau=LAG / fs)[0] * fs) == TRUE
    c += plain_cc(m2, m1, LAG) == TRUE
print(f"echo      : GCC-PHAT {p}/20  plain {c}/20")

# 実験2：各マイクに独立な白色雑音を足した場合
for snr in (10, 0, -5, -10):
    p = c = 0
    for s in range(50):
        rng = np.random.default_rng(100 + s)
        src = source(rng)
        sd = np.sqrt(np.var(src) / 10 ** (snr / 10))
        m1 = src + rng.standard_normal(N) * sd
        m2 = np.roll(src, TRUE) + rng.standard_normal(N) * sd
        p += round(gcc_phat(m2, m1, fs, max_tau=LAG / fs)[0] * fs) == TRUE
        c += plain_cc(m2, m1, LAG) == TRUE
    print(f"SNR {snr:>3} dB: GCC-PHAT {p}/50  plain {c}/50")

# 遅れ → 到来角（2マイク・間隔10cm・音速343m/s）
d, v = 0.10, 343.0
tau = TRUE / fs
print(f"tau={tau * 1e3:.4f} ms -> {np.degrees(np.arcsin(tau * v / d)):.1f} deg")

実行結果は次のとおりです。

echo      : GCC-PHAT 20/20  plain 14/20
SNR  10 dB: GCC-PHAT 50/50  plain 50/50
SNR   0 dB: GCC-PHAT 20/50  plain 24/50
SNR  -5 dB: GCC-PHAT 11/50  plain 11/50
SNR -10 dB: GCC-PHAT 5/50  plain 9/50
tau=0.1875 ms -> 40.0 deg

最終行は、3サンプル（0.1875ms）の遅れが、間隔10cmのマイクペアでは約40.0°の到来角に当たることを示しています。

補間倍率 interp と探索範囲 max_tau の決め方

interp=16 は、逆FFTの長さを16倍にしてサンプル間を1/16刻みで評価する設定です。16kHz・10cm間隔では、遅れが0から1サンプルに変わると角度は約12.4°変わります。角度刻みは正面から離れるほど大きくなるため、補間なしでは角度が不均等な階段状になります。max_tau にはマイク間隔を音速で割った値を渡してください。物理的にあり得ない遅れを探索範囲から外すだけで、反射音による誤検出が減ります。

フィルタ設計やリサンプリングを組み合わせる場合は、SciPyとは？読み方・使い方・NumPyとの違いをPythonで解説の scipy.signal が使えます。

MATLABのgccphat関数の構文と戻り値

MATLABでは、Phased Array System Toolboxの gccphat 関数（R2015bで導入）が同じ処理を提供しています。tau = gccphat(sig,refsig,fs) で基準信号に対する遅れを秒で返し、[tau,R,lag] = gccphat(___) とすれば相関値とラグ時間も受け取れます。MATLAB本体とツールボックスの関係はMATLAB/Simulinkとは？何ができる・違い・料金を初心者向けに解説を参照してください。

反射音と雑音でPHATの効き方を比べた検証結果

上のコードでは、低域に偏った信号（白色雑音の64サンプル移動平均）を使い、真の遅れ3サンプルを当てられた回数を数えています。探索範囲は比較実験用に両方式とも±40サンプルとし、正解は推定遅れを整数サンプルに丸めて3になる場合としました。10cm間隔の実マイクでは、この範囲を約±4.66サンプル以内に制限します。遅延生成には配列末尾を先頭へ回すnp.rollを使っており、実室内の伝搬を再現した実験ではありません。乱数の種を変えた20回または50回の結果は、この合成条件での傾向を示すものです。

単一反射音を加えた条件でのPHATの正解数

直接音の20サンプル後に振幅0.8の反射音を足すと、GCC-PHATは20回中20回正解し、重みなしの相互相関は14回にとどまりました。低域に偏った信号では重みなしの相関ピークが幅広くなり、反射音のピークと重なって最大値の位置がずれるためです。PHATで白色化するとピークが細くなり、直接音と反射音が分離します。

白色雑音を加えた低SNR条件での正解数比較

各マイクに独立な白色雑音を足すと、SNR 10dBでは両方式とも50回全問正解でした。SNR 0dBではGCC-PHATが20回、重みなしが24回で、-10dBでは5回対9回と逆転しました。この信号は高域にほとんどエネルギーが無いため、高域の成分はほぼ雑音だけです。PHATは雑音しかない帯域まで振幅1に持ち上げて、信号のある帯域と同じ重みで足し合わせてしまいます。

今回の合成信号では、PHATは追加した反射音による誤推定を減らしましたが、独立な白色雑音下で重みなしの相互相関を上回りませんでした。定常雑音が大きい環境では、前処理で雑音を抑えるか、重みの設計を見直してください。

SRP-PHATとGCC-PHATの違い

SRP-PHAT（Steered Response Power with Phase Transform）は、候補となる位置や方向を空間に並べ、それぞれについて全マイクペアのGCC-PHAT値を足し合わせ、合計が最大になる候補を音源位置とする方式です。DiBiase、Silverman、Brandsteinが2001年の書籍『Microphone Arrays』（Springer）の章「Robust Localization in Reverberant Rooms」で体系的に示しました。

ペアごとの推定をあとで幾何計算でまとめる方式は軽いものの、1ペアの誤ピークがそのまま結果を崩します。SRP-PHATは誤ピークに強い代わりに、候補点の数だけ計算が増えます。マイクが2本ならGCC-PHAT、4本以上で3次元の位置まで必要ならSRP-PHAT、という分け方が出発点になります。

GCC-PHATの実装事例：GoogleのSpeechCompass

Googleの研究チームがCHI 2025で発表したSpeechCompassは、スマートフォン向けの字幕表示で、話者の方向を矢印や色分けで示すシステムです。論文によると、4本のデジタルマイク（STMicroelectronics MP34DT01-M）と110MHzのArm Cortex-M33マイコン（STM32L55）を使い、GCC-PHATの変種をArmのCMSISライブラリで実装しています。

1ペアの遅れ計算に2.9msかかり、6ペア分のデータを揃えるのに17.4msとされています。推定値は直近600サンプルに対するガウスカーネル密度推定でまとめ、通常の会話音量（60〜65dB）での誤差は11.1〜22.1°でした。17.4msは全ペアの計算時間であり、安定した方向推定までの時間とは異なります。同論文では、音声の方向推定に平均約263ms、音声の特性により約70〜500msを要したと報告しています。

GCC-PHATを選ばないほうがよい場面

• 複数の話者が同時に話す：GCC-PHATは各フレームで最も強い1つのピークを採る設計です。同時発話では相関関数に複数のピークが立ち、どれが目的の話者かを区別できません。複数話者の方向推定には、MUSICや複数ピークを探索するSRP-PHATの構成を検討してください。話者数の推定や音声波形の分離には、別途そのための処理が必要です。
• 定常雑音が支配的でSNRが低い：前章の低域に偏った合成信号と独立な白色雑音の検証では、SNR 0、-5、-10dBで重みなしの相互相関と同等以下の正解数でした。
• 音源の種類ごとに位置を出したい：音響イベントの検出と定位を同時に行うSELDでは、GCC-PHATを単独の推定器にせず、ニューラルネットワークへ入力する特徴量として使う研究があります。

よくある質問

GCC-PHATのPHATは何の略ですか？

Phase Transform（位相変換）の略です。相互スペクトルを振幅で割り、位相の情報だけを残す重み付けを指します。

GCC-PHATとSRP-PHATの違いは何ですか？

GCC-PHATは2本のマイクの時間差を1つ求める手法です。SRP-PHATは空間上の候補点ごとに全ペアのGCC-PHAT値を合算し、最大になる位置を探します。

PythonでGCC-PHATを使うにはどうすればよいですか？

NumPyの np.fft.rfft と np.fft.irfft を使えば、本記事の関数のように十数行で実装できます。探索範囲はマイク間隔÷音速に制限してください。

正規化相互相関とGCC-PHATは同じものですか？

同じではありません。正規化相互相関は全体のエネルギーで1回だけ割り、GCC-PHATは周波数ごとに振幅で割ります。

GCC-PHATは雑音に強いのですか？

強いのは主に残響（反射音）に対してです。本記事の検証では、SNR 0dB以下の白色雑音下で重みなしの相互相関と同等以下の正解率でした。

関連記事

• NumPyとは？読み方・インストール・基本的な使い方をPythonで解説
• SciPyとは？読み方・使い方・NumPyとの違いをPythonで解説
• テンプレートマッチングとは？matchTemplateの仕組みと類似度指標の選び方を解説
• カルマンフィルタとは？仕組みと使い方をわかりやすく解説【線形の基礎から応用まで】
• MATLAB/Simulinkとは？何ができる・違い・料金を初心者向けに解説