GCC-PHATとは？PHAT重みの意味とTDOA推定の手順・Python実装と残響室での精度

GCC-PHAT（Generalized Cross-Correlation with Phase Transform）は、2本のマイクに届いた音の時間差（TDOA：Time Difference of Arrival）を求めるアルゴリズムです。相互スペクトルの振幅を周波数ごとに1へそろえ、位相だけを残してから時間領域へ戻すため、相関のピークが鋭くなります。この記事では、PHATという重みの意味、計算の手順、NumPyでの実装を説明したうえで、Pythonの室内音響シミュレーションで残響時間を変えながら精度を測った結果を示します。PHATが効く残響の範囲、効かなくなる音源の種類、重み指数βを下げたときの変化も、実行した数値で確認できます。

まとめ：GCC-PHATの要点

• GCC-PHATは、一般化相互相関（GCC）の重み関数に「相互スペクトルの振幅の逆数」を使う方式です。GCCの枠組みはKnappとCarterが1976年に示しました。
• PHATは位相変換（Phase Transform）の略です。振幅を捨てて位相だけで相関を取るため、反射音のピークと直接音のピークが分かれやすくなります。
• 手順は「フレーム切り出し→FFT→相互スペクトル→振幅で割る→逆FFT→ピーク探索」の6段で、NumPyなら関数1つ・十数行で書けます。
• pyroomacoustics 0.10.1で6×5×3mの部屋を再現し、残響時間を目標0.6秒（実測0.76秒）、SNR 20dBにした条件では、到来角の誤差5°以内に収まった割合が重みなしの相互相関で47.5%、GCC-PHATで66.7%でした。8フレーム分の相関関数を足してからピークを取ると73.3%に上がります。
• 今回のシミュレーションでは、残響が短い（実測0.33秒）か無い部屋で重みなしとの差はほぼありませんでした。基本周波数120〜220Hzの調波音のような狭帯域の音では、GCC-PHATの正解率は重みなしを下回りました。
• 多数のマイクで空間を走査する場合は、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」（IEEE Transactions on Acoustics, Speech, and Signal Processing 24巻4号、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)| です。相互スペクトルの大きさが0でない周波数では、割った結果の大きさは1になります。両信号が時間差τ0だけを持つ理想条件なら、残るのは位相 e^(-j2πfτ0) だけになります。これを逆変換した結果は τ0 の位置に立つ鋭いインパルスになります。

重みなしの相互相関は、エネルギーの大きい周波数帯に結果が引っ張られます。音声のエネルギーは低い周波数帯に偏るため、相関のピークがなだらかに広がり、近くに反射音のピークがあると区別しにくくなります。PHATはスペクトルを平らにする白色化の一種で、電子情報通信学会の知識ベース（2群6編3章）では、クロススペクトルの振幅で白色化する同じ処理を「白色化相互相関（Cross-power Spectrum Phase Analysis：CSP）法」と呼んでいます。音響信号処理で使う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ペアです。複数ペアの時間差から位置を出す計算は、TDOA測位（双曲線測位）と呼ばれます。各ペアのτは音源が乗る双曲線（3次元では双曲面）を1つ決めるので、複数ペアの交点を最小二乗で求めます。

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

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

動く音源を追う場合、フレームごとの推定値はばらつきます。実用システムでは推定値をそのまま使わず、ヒストグラムやカルマンフィルタとは？仕組みを1次元の例でわかりやすく解説【Python実装とQ・Rの決め方】のような追跡処理で平滑化します。止まっている音源なら、後述の検証のように複数フレームの相関関数を足してからピークを取る方法も効きます。

PythonによるGCC-PHAT実装

NumPyで書くgcc_phat関数

次の関数は、2026年10月時点でGitHubのスターが200を超えるxiongyihui/tdoaの gcc_phat.py（Apache License 2.0）と同じ構成で、重み指数 beta と0除算対策を加え、ピーク探索を相関値の絶対値で行うように変えたものです。依存はNumPyとは？読み方・インストール・基本的な使い方をPythonで解説だけです。

import numpy as np

# 関数の構成は xiongyihui/tdoa の gcc_phat.py（Apache License 2.0）を参考にし、
# 重み指数 beta・0除算対策・探索範囲の切り出しを加えた
def gcc_phat(sig, refsig, fs, max_tau, interp=16, beta=1.0):
    """refsig を基準にした sig の遅れ（秒）と相関関数を返す。beta=0で重みなし、1でPHAT"""
    n = len(sig) + len(refsig)
    R = np.fft.rfft(sig, n=n) * np.conj(np.fft.rfft(refsig, n=n))  # 相互スペクトル
    R = R / (np.abs(R) ** beta + 1e-12)                             # PHAT-β重み
    cc = np.fft.irfft(R, n=interp * n)                              # 零詰めで補間しながら時間領域へ
    m = min(int(interp * fs * max_tau), int(interp * n / 2))
    cc = np.concatenate((cc[len(cc) - m:], cc[:m + 1]))             # ±max_tau の範囲だけ残す
    if not np.any(cc):
        return np.nan, cc
    shift = np.argmax(np.abs(cc)) - m
    return shift / (interp * fs), cc

使い方は tau, cc = gcc_phat(mic2, mic1, 16000, max_tau=0.10 / 343) のように、基準にするマイクの信号を2番目に渡します。戻り値の tau が正なら、1番目の信号が後から届いています。beta=1.0 が通常のGCC-PHAT、beta=0.0 が重みなしの相互相関です。間の値については後述のPHAT-βの章で扱います。入力が無音で相関値がすべて0のときは、ピーク位置が決まらないため tau に nan を返します。

補間倍率 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(___) とすれば相関値とラグ時間も受け取れます。参照信号を省いた gccphat(sig,fs) はM本のチャネルの全組み合わせを計算し、tauは1×M²の行になります。サンプリング周波数を省くと既定値1Hzで計算されるため、秒単位の遅れが欲しいときは fs を必ず渡してください。MATLAB本体とツールボックスの関係はMATLAB/Simulinkとは？何ができる・違い・料金を初心者向けに解説を参照してください。

残響室シミュレーションで測ったGCC-PHATの精度

PHATがどの程度の残響まで効くのかを確かめるため、Pythonの室内音響シミュレーションライブラリpyroomacoustics（0.10.1、MITライセンス）の鏡像法で部屋の反射を再現し、重みなし（β=0）、PHAT-β（β=0.7）、GCC-PHAT（β=1）の3つを比べました。

• 部屋：6×5×3mの直方体。反射なし（max_order=0で直接音のみ）と、inverse_sabine に目標残響時間0.3、0.6、0.9秒を渡して壁の吸音率を決めた3条件
• マイク：間隔10cmの2本を部屋の中央、高さ1.2mに配置。音源は1.5m先、到来角を-60°〜60°から一様に抽選
• 音源：パワーが周波数に反比例するピンクノイズ（広帯域）と、基本周波数120〜220Hz・第5倍音までの調波音（狭帯域）
• 雑音：各マイクに独立な白色雑音を加算。SNR 20dB・0dBは、2000サンプル目以降の録音全体（残響の減衰部を含む）の分散を基準にした。評価に使う8フレームだけで測ると実効SNRはこれより高く、1配置で測った例では目標0.6秒で約3.9dB、0.9秒で約5.0dB上振れする
• 評価：16kHz、1フレーム1024サンプル（64ms）。1条件30配置×8フレームで、推定角の誤差が5°以内に入った割合

精度評価とは別の代表配置（音源座標［3.0, 4.0, 1.2］m、マイク座標［3.0, 2.5, 1.2］m）でインパルス応答から測った残響時間は、目標0.3秒に対して0.33秒、0.6秒に対して0.76秒、0.9秒に対して1.14秒でした。本記事の「実測」はこの代表配置での参考値であり、精度評価の全配置・両マイクで測った値ではありません。以下の表では実測値を併記します。

pyroomacoustics 0.10.1による検証コード

前章の gcc_phat 関数の後ろに次のコードを続けて実行します。pip install pyroomacoustics==0.10.1 で入ります。Python 3.9.6・NumPy 2.0.2・pyroomacoustics 0.10.1で実行し、所要時間は約7〜9分でした。

import pyroomacoustics as pra

FS, C, D = 16000, 343.0, 0.10        # サンプリング周波数・音速・マイク間隔
MAX_TAU = D / C                      # 物理的にあり得る最大の遅れ（約0.29ms）
FRAME, NFRAME, TRIALS = 1024, 8, 30  # 64msの窓を8枚、条件ごとに30配置
BETAS = (0.0, 0.7, 1.0)

def pink(n, rng):
    """パワーが周波数に反比例するピンクノイズ（広帯域の例）"""
    X = np.fft.rfft(rng.standard_normal(n))
    f = np.arange(len(X), dtype=float)
    f[0] = 1.0
    return np.fft.irfft(X / np.sqrt(f), n)

def harmonic(n, rng):
    """基本周波数120〜220Hz・第5倍音までの調波音（狭帯域の例）"""
    t = np.arange(n) / FS
    f0 = rng.uniform(120, 220)
    return sum(np.sin(2 * np.pi * k * f0 * t + rng.uniform(0, 2 * np.pi)) / k for k in range(1, 6))

def simulate(rt60, signal, rng):
    """6×5×3mの部屋で、2本のマイクの1.5m先に音源を置いて録音を合成する"""
    dim = [6.0, 5.0, 3.0]
    if rt60 > 0:
        e_abs, max_order = pra.inverse_sabine(rt60, dim)
        room = pra.ShoeBox(dim, fs=FS, materials=pra.Material(e_abs), max_order=max_order)
    else:
        room = pra.ShoeBox(dim, fs=FS, max_order=0)  # 反射なし（直接音のみ）
    center = np.array([3.0, 2.5, 1.2])
    ang = rng.uniform(-60, 60)                       # 真の到来角（度）
    rad = np.radians(ang)
    room.add_source(list(center + 1.5 * np.array([np.sin(rad), np.cos(rad), 0])), signal=signal)
    room.add_microphone_array(np.c_[center - [D / 2, 0, 0], center + [D / 2, 0, 0]])
    room.simulate()
    return room.mic_array.signals, ang

def to_deg(tau):
    return np.degrees(np.arcsin(np.clip(tau * C / D, -1, 1)))

for name, gen, rt60s in (("pink", pink, (0.0, 0.3, 0.6, 0.9)), ("harmonic", harmonic, (0.3, 0.6))):
    for rt60 in rt60s:
        for snr in (20, 0):
            rng = np.random.default_rng(int(rt60 * 10) * 100 + snr)
            frame_ok = {b: [] for b in BETAS}
            sum_ok = {b: [] for b in BETAS}
            for _ in range(TRIALS):
                s, ang = simulate(rt60, gen(FS, rng), rng)
                noise_sd = np.sqrt(np.var(s[0][2000:]) / 10 ** (snr / 10))
                s = s + rng.standard_normal(s.shape) * noise_sd  # 各マイクに独立な白色雑音
                for b in BETAS:
                    total = 0
                    for st in range(2000, 2000 + NFRAME * FRAME, FRAME):
                        tau, cc = gcc_phat(s[0][st:st + FRAME], s[1][st:st + FRAME], FS, MAX_TAU, beta=b)
                        frame_ok[b].append(abs(to_deg(tau) - ang) <= 5)
                        total = total + cc               # 8フレームの相関関数を足し合わせる
                    m = int((len(total) - 1) / 2)
                    sum_ok[b].append(abs(to_deg((np.argmax(np.abs(total)) - m) / (16 * FS)) - ang) <= 5)
            print(f"{name:8s} RT60={rt60:.1f}s SNR={snr:>2}dB", " ".join(
                f"| beta={b}: {np.mean(frame_ok[b]) * 100:5.1f}% sum8 {np.mean(sum_ok[b]) * 100:5.1f}%" for b in BETAS), flush=True)

# inverse_sabine に渡した目標値と、生成したインパルス応答から測った残響時間
for rt60 in (0.3, 0.6, 0.9):
    e_abs, max_order = pra.inverse_sabine(rt60, [6.0, 5.0, 3.0])
    room = pra.ShoeBox([6.0, 5.0, 3.0], fs=FS, materials=pra.Material(e_abs), max_order=max_order)
    room.add_source([3.0, 4.0, 1.2])
    room.add_microphone([3.0, 2.5, 1.2])
    room.compute_rir()
    print(f"target RT60={rt60:.1f}s -> measured {room.measure_rt60()[0][0]:.2f}s")

実行結果は次のとおりです。各行の左の値が1フレームごとの正解率、sum8 が8フレーム分の相関関数を足してから推定した場合の正解率です。1フレームの値は30配置×8フレーム＝240回、sum8は30配置での割合で、1配置が3.3ポイントに当たります。

pink     RT60=0.0s SNR=20dB | beta=0.0: 100.0% sum8 100.0% | beta=0.7: 100.0% sum8 100.0% | beta=1.0: 100.0% sum8 100.0%
pink     RT60=0.0s SNR= 0dB | beta=0.0:  95.8% sum8 100.0% | beta=0.7:  95.4% sum8 100.0% | beta=1.0:  93.8% sum8 100.0%
pink     RT60=0.3s SNR=20dB | beta=0.0:  98.8% sum8 100.0% | beta=0.7: 100.0% sum8 100.0% | beta=1.0: 100.0% sum8 100.0%
pink     RT60=0.3s SNR= 0dB | beta=0.0:  80.8% sum8  96.7% | beta=0.7:  80.0% sum8  96.7% | beta=1.0:  78.8% sum8 100.0%
pink     RT60=0.6s SNR=20dB | beta=0.0:  47.5% sum8  46.7% | beta=0.7:  65.8% sum8  70.0% | beta=1.0:  66.7% sum8  73.3%
pink     RT60=0.6s SNR= 0dB | beta=0.0:  49.6% sum8  46.7% | beta=0.7:  53.3% sum8  70.0% | beta=1.0:  50.0% sum8  73.3%
pink     RT60=0.9s SNR=20dB | beta=0.0:  40.4% sum8  36.7% | beta=0.7:  53.8% sum8  56.7% | beta=1.0:  54.2% sum8  63.3%
pink     RT60=0.9s SNR= 0dB | beta=0.0:  29.6% sum8  26.7% | beta=0.7:  37.1% sum8  33.3% | beta=1.0:  33.8% sum8  36.7%
harmonic RT60=0.3s SNR=20dB | beta=0.0:  23.3% sum8  23.3% | beta=0.7:   7.5% sum8   3.3% | beta=1.0:   6.7% sum8   0.0%
harmonic RT60=0.3s SNR= 0dB | beta=0.0:  16.2% sum8  20.0% | beta=0.7:  10.8% sum8  10.0% | beta=1.0:   9.6% sum8  10.0%
harmonic RT60=0.6s SNR=20dB | beta=0.0:  17.1% sum8  20.0% | beta=0.7:  11.7% sum8  10.0% | beta=1.0:  10.4% sum8  10.0%
harmonic RT60=0.6s SNR= 0dB | beta=0.0:  13.8% sum8  23.3% | beta=0.7:  10.4% sum8  13.3% | beta=1.0:   9.2% sum8  13.3%
target RT60=0.3s -> measured 0.33s
target RT60=0.6s -> measured 0.76s
target RT60=0.9s -> measured 1.14s

残響時間ごとの正解率（広帯域のピンクノイズ）

ピンクノイズ音源の1フレームごとの正解率を、残響時間とSNRで並べると次のとおりです。

差が開いたのは残響が長く、雑音が小さい条件です。実測0.76秒・SNR 20dBでは、GCC-PHATが重みなしより19.2ポイント高くなりました。今回のピンクノイズ音源では、反射なしと目標残響時間0.3秒の条件で3方式とも78.8%以上となり、1フレーム評価の方式間の差は2ポイント以内でした。この条件ではPHATによる改善は小さいものの、音声や別のマイク配置でも同じ結果になるとは限りません。反対に、SNR 0dBでは残響0.76秒でも差が0.4ポイントまで縮みます。雑音が支配的な帯域もPHATが同じ重みで扱うことは、性能低下の一因と考えられます。ただし、このコードではSNRごとに乱数シードを変えており、音源波形と配置も異なるため、条件間の差を雑音の増加だけに帰することはできません。

8フレームの相関関数を合算する効果

今回の広帯域ピンクノイズに対するGCC-PHATでは、相関関数を8フレーム分足してからピークを取ると、1フレーム評価以上の正解率になりました。ただし、音源が静止していても、方式や音源によっては合算後に正解率が下がります。実測0.76秒・SNR 0dBでは、GCC-PHATの1フレームの50.0%が、8フレーム（512ms）の合算で73.3%になりました。重みなしは同じ条件で1フレームの49.6%から合算後46.7%へ下がり、合算しても伸びません。重みなしの相関では、反射音による偽のピークが毎フレーム同じ位置に出るため、足しても偽のピークが一緒に育つと考えられます。

合算できるのは、合算する時間内に音源の方向がほぼ変わらない場合に限られます。移動する話者を追うなら合算を短くし、追跡処理側で平滑化してください。

狭帯域の調波音でPHATが負ける条件

基本周波数120〜220Hz、第5倍音まで（最大で約1.1kHz）の調波音では、どの方式も正解率が23.3%以下にとどまりました。残響0.33秒・SNR 20dBでは、重みなしが23.3%、β=0.7が7.5%、GCC-PHATが6.7%で、8フレーム合算のGCC-PHATは0.0%です。誤差5°以内の幅10°は抽選範囲120°の約8.3%に当たり、GCC-PHATの6.7〜10.4%はでたらめに角度を選ぶのと変わらない水準です。この音は数本の倍音以外の帯域にエネルギーがほとんど無く、PHATはそこにある雑音と残響の成分まで倍音と同じ重みで足し合わせます。

200Hzの波長は約1.7mあり、間隔10cmのマイクでは位相差がわずかしか生じません。モーター音やハムのような狭帯域の音源を定位したいなら、PHATの重みを調整するより先に、マイク間隔を広げるか、狭帯域向けの方式を選ぶべきです。

PHAT-β：振幅情報を一部残す重み指数

PHAT-βは、相互スペクトルの振幅情報をどこまで残すかを1つのパラメータβで調整する変形です。Donohue、Hannemann、Dietzの論文「Performance of phase transform for detecting sound sources with microphone arrays in reverberant and noisy environments」（Signal Processing 87巻7号、2007年、1677〜1691ページ）が提案しました。8素子のアレイで遅延和の前にPHATを掛けるモンテカルロシミュレーションを評価した要旨では、通常のPHATは広帯域の信号で検出性能を大きく改善する一方、狭帯域の目標では性能を落としやすく、PHAT-βは両方で改善を得られるとしています。実音声での結果から、一般の信号を扱うアレイ用途ではβを0.5〜0.7にすることを提案しています。

本記事の関数では、相互スペクトルを振幅のβ乗で割る形で実装しました（β=0で重みなし、β=1で通常のPHAT）。今回のシミュレーションでβ=0.7は、ピンクノイズではGCC-PHATとほぼ同じ正解率で、実測0.76秒・SNR 0dBの1フレーム評価では53.3%とGCC-PHATの50.0%を上回りました。調波音では重みなしとGCC-PHATの中間で、重みなしを上回る条件はありませんでした。

実務での判断は次のとおりです。会議室の話し声のように音源の帯域が広く、残響が主な敵ならβ=1のGCC-PHATで始めます。空調などの定常雑音が大きい現場や、音源の帯域が分からない汎用の装置では、β=0.5〜0.7を候補にして実環境の録音で比べてください。βだけで狭帯域の音源の定位を救うことは、今回の条件ではできませんでした。

SRP-PHATとGCC-PHATの違い

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

ペアごとの推定をあとで幾何計算でまとめる方式は軽いものの、1ペアの誤ピークがそのまま結果を崩します。SRP-PHATは誤ピークに強い代わりに、候補点の数だけ計算が増えます。マイクが2本ならGCC-PHAT、4本以上で3次元の位置まで必要ならSRP-PHAT、という分け方が出発点になります。2024年に公開されたGrinsteinらのチュートリアルレビュー（arXiv:2405.02991）は、SRP-PHATに関する200本超の論文を計算量の削減、雑音・残響への頑健化、複数音源への拡張の観点で整理し、モジュール化した実装X-SRPをPythonで公開しています。

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°でした。音が小さくなるほど誤差は大きく増え、論文はその理由を、GCC-PHATが最も大きな音を拾うため、音源が環境雑音と同じかそれ以下の大きさだと時間差を推定しにくいと説明しています。17.4msは全ペアの計算時間であり、安定した方向推定までの時間とは異なります。同論文では、音声の方向推定が安定するまで平均約263ms、音声の性質により約70〜500msかかったと報告しています。

同じ論文は、市販のPixel 6でも試しています。この機種はマイクを3本持ちますが、2本の間隔が15mmと近すぎて、使えるペアが2組に限られるため、定位は180°の範囲に限られました。平均誤差は音声で16.0°、環境音で23.0°でした。

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

• 複数の話者が同時に話す：GCC-PHATは各フレームで最も強い1つのピークを採る設計です。同時発話では相関関数に複数のピークが立ち、どれが目的の話者かを区別できません。複数話者の方向推定には、MUSICや複数ピークを探索するSRP-PHATの構成を検討してください。
• 音源が狭帯域：前述の調波音の検証では、GCC-PHATの正解率は重みなしの相互相関を下回りました。
• 長い残響下で単一ペアでは要求精度を満たせない：実測1.14秒・SNR 20dBでは、GCC-PHATでも1フレームの正解率は54.2%でした。マイクを音源に近づける、吸音するなどの収音環境の改善に加え、マイクを追加して複数ペアを使うSRP-PHATへ切り替える方法も検討します。
• 定常雑音が支配的でSNRが低い：SNR 0dBでは、残響0.76秒でもGCC-PHATと重みなしの差は0.4ポイントでした。雑音抑圧の前処理かβの調整を先に検討します。

音源定位を組み込んだ装置や業務システムを、既製品ではなく要件に合わせて作る場合は、フルスクラッチ開発の範囲で、マイク配置の設計から検証まで相談できます。

よくある質問

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 を使えば、本記事の関数のように十数行で実装できます。探索範囲はマイク間隔÷音速に制限してください。

TDOAで音源の位置を推定するには何本のマイクが必要ですか？

平面内の方向だけなら2本で求まりますが、前後の区別がつきません。360°の方位には同一直線上にない3本以上、3次元の位置には同一平面上にない4本以上が目安です。

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

強いのは主に残響（反射音）に対してです。本記事の検証では、SNR 0dBの白色雑音を加えると、残響0.76秒でも重みなしの相互相関との差は0.4ポイントでした。

関連記事

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