基于 GCC-PHAT 的 TDOA 估计方法

最新推荐文章于 2026-05-28 16:00:59 发布

原创最新推荐文章于 2026-05-28 16:00:59 发布 · 1k 阅读

23 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

声源定位算法

在前文中已经系统介绍了 CBF、MVDR、MUSIC、CLEAN-SC 等空间谱类方法。这些方法本质上基于阵列协方差矩阵与空间导向矢量进行建模与反演。而在阵列信号处理中，还有一类同样重要但思想不同的方法——基于到达时间差（TDOA, Time Difference of Arrival）的定位方法。

TDOA 是阵列定位中最基本的观测量之一。在远场假设下，两麦之间的时间差即可直接映射为方位角；在近场或多阵元情况下，多个 TDOA 可构成双曲线交汇定位问题。GCC-PHAT 是当前工程中最常用、最稳健的 TDOA 估计方法之一。

本文系统给出其数学推导、工程实现逻辑及物理约束条件。

一、TDOA 问题的数学模型

设空间中存在单一声源信号 s(t)，两只麦克风接收信号为：

$x_1(t) = s(t) + n_1(t)$

$x_2(t) = s(t-\tau) + n_2(t)$

其中：

τ 为两麦之间的传播时延
n1(t),n2(t)为噪声
假设直达声占主导

离散形式为：

$x_2[n] = x_1[n - \ell_0]$

其中：

$\tau = \frac{\ell_0}{f_s}$

二、互相关法的基本思想

TDOA 的最直接估计方式是互相关函数：

$R_{12}(\ell) = \sum_{n} x_1[n] x_2[n+\ell]$

估计延迟：

$\hat{\ell} = \arg\max_{\ell} R_{12}(\ell)$

$\hat{\tau} = \frac{\hat{\ell}}{f_s}$

问题

直接互相关在实际工程中存在严重缺陷：

低频或强能量频段主导相关结果
混响环境下峰值变宽
窄带信号导致多峰模糊

因此需要频域加权改进。

三、广义互相关（GCC）

互相关可以在频域中表示为：

$G_{12}(k) = X_1(k) X_2^*(k)$

$R_{12}(\ell) = \text{IFFT}\{G_{12}(k)\}$

广义互相关形式为：

$R^{GCC}_{12}(\ell) = \text{IFFT}\{\Psi(k) G_{12}(k)\}$

其中 Ψ(k)为频域权重函数。

四、PHAT 权重

PHAT 权重定义为：

$\Psi_{PHAT}(k) = \frac{1}{|G_{12}(k)| + \epsilon}$

因此：

$R^{PHAT}_{12}(\ell) = \text{IFFT}\left\{ \frac{G_{12}(k)} {|G_{12}(k)| + \epsilon} \right\}$

PHAT 的核心思想

将互谱幅值归一化，仅保留相位信息。

其物理意义是：

延迟信息体现在相位差
幅度结构在混响环境中高度失真
相位一致性比能量一致性更稳定

因此 PHAT 本质上是一种“相位对齐增强”。

五、物理约束条件

两麦间距为 d，声速为 c：

$|\tau| \le \frac{d}{c}$

远场条件下：

$\tau = \frac{d \sin\theta}{c}$

因此：

$\theta = \arcsin\left(\frac{c\tau}{d}\right)$

六、与空间谱类方法的关系

方法	核心思想	是否依赖协方差矩阵
CBF	波束形成	是
MVDR	最小方差约束	是
MUSIC	子空间分解	是
CLEAN-SC	迭代去卷积	是
GCC-PHAT	相位对齐求延迟	否

GCC-PHAT 与上述方法的区别在于：

它不构造空间谱
不做特征值分解
不显式建模信号子空间

而是直接估计几何延迟。

从工程角度看：

GCC-PHAT 更轻量
对实时系统友好
适合低通道数阵列

但在多源分辨能力方面弱于 MUSIC 等子空间方法。

七、工程实现关键点

1）分帧处理

避免长时间非平稳信号破坏相关结构。

2）频段限制

带通滤波能显著提高稳定性。

3）亚采样插值

抛物线拟合或频域零填充提高估计精度。

4）多帧稳健统计

可采用中值滤波或 RANSAC 抑制跳变。

八、仿真实现

代码：

%% ==============================================================
%  GCC-PHAT TDOA 仿真对比
%  1) xcorr
%  2) MATLAB 内置 gccphat
%  3) 手动 GCC-PHAT
% ==============================================================

clear; clc; close all;

%% =========================
% 1) 参数设置
% =========================
fs   = 48000;          % 采样率
T    = 0.25;           % 信号长度（秒）
N    = round(T*fs);

trueDelay = 37;        % 真实延迟（样本点）
SNRdB = 15;            % 加噪信噪比

f_full_lo  = 20;       % 生成信号的全频带
f_full_hi  = 20000;

f_focus_lo = 300;      % 关注频段
f_focus_hi = 8000;

epsPHAT = 1e-12;

%% =========================
% 2) 生成 20Hz~20kHz 带限宽带信号
% =========================
x = randn(N,1);

nfftGen = 2^nextpow2(N);
X = fft(x,nfftGen);

f = (0:nfftGen-1)'*(fs/nfftGen);
mask = false(nfftGen,1);

mask(f>=f_full_lo & f<=f_full_hi) = true;
mask(f>=fs-f_full_hi & f<=fs-f_full_lo) = true;

X = X .* mask;
x1 = real(ifft(X,nfftGen));
x1 = x1(1:N);
x1 = x1 / rms(x1);

%% =========================
% 3) 构造第二通道 + 加噪
% =========================
x2 = [zeros(trueDelay,1); x1(1:end-trueDelay)];

sigPow = mean(x1.^2);
noiPow = sigPow / (10^(SNRdB/10));

x1n = x1 + sqrt(noiPow)*randn(size(x1));
x2n = x2 + sqrt(noiPow)*randn(size(x2));

%% =========================
% 4) 300–8000 Hz 带通
% =========================
bp = designfilt('bandpassiir', ...
    'FilterOrder',8, ...
    'HalfPowerFrequency1',f_focus_lo, ...
    'HalfPowerFrequency2',f_focus_hi, ...
    'SampleRate',fs);

x1bp = filtfilt(bp,x1n);
x2bp = filtfilt(bp,x2n);

%% =========================
% 5) 方法A：xcorr
% =========================
[r_xcorr,lags_xcorr] = xcorr(x2bp,x1bp);
[~,idxA] = max(r_xcorr);
lagA = lags_xcorr(idxA);

%% =========================
% 6) 方法B：内置 gccphat
% =========================
[tau_est,r_builtin,lags] = gccphat(x2bp,x1bp,fs);

lagB = round(tau_est * fs);   % 秒 -> 样本
lag_builtin = lags * fs;      

%% =========================
% 7) 方法C：手动 GCC-PHAT（频域门控）
% =========================
nfft = 2^nextpow2(2*N-1);

X1 = fft(x1n,nfft);
X2 = fft(x2n,nfft);

G12 = X2 .* conj(X1);

fC = (0:nfft-1)'*(fs/nfft);
maskFocus = false(nfft,1);

maskFocus(fC>=f_focus_lo & fC<=f_focus_hi) = true;
maskFocus(fC>=fs-f_focus_hi & fC<=fs-f_focus_lo) = true;

G12 = G12 .* maskFocus;          % 频域门控
G12 = G12 ./ (abs(G12)+epsPHAT); % PHAT

r_temp = ifft(G12,nfft,'symmetric');
r_temp = fftshift(r_temp);

mid = floor(nfft/2)+1;
idxRange = (mid-(N-1)):(mid+(N-1));

r_manual = r_temp(idxRange);
lags_manual = (-(N-1):(N-1))';

[~,idxC] = max(r_manual);
lagC = lags_manual(idxC);

%% =========================
% 8) 打印结果
% =========================
fprintf('\n========== 延迟估计结果 ==========\n');
fprintf('True delay (samples)      : %d\n', trueDelay);
fprintf('xcorr (bandpass)          : %d\n', lagA);
fprintf('builtin gccphat (bandpass): %d\n', lagB);
fprintf('manual GCC-PHAT (mask)    : %d\n', lagC);

%% =========================
% 9) 画完整相关函数
% =========================
figure;
plot(lags_xcorr,r_xcorr./max(abs(r_xcorr)),'LineWidth',1); hold on;
plot(lag_builtin,r_builtin./max(abs(r_builtin)),'LineWidth',1);
plot(lags_manual,r_manual./max(abs(r_manual)),'LineWidth',1);

grid on;
xlabel('Lag (samples)');
ylabel('Normalized Correlation');
title('相关函数对比');
legend('xcorr(带通)','gccphat(带通)','手动GCC-PHAT(频域门控)','Location','best');

%% =========================
% 10) 峰值附近放大
% =========================
win = 50;
lo = trueDelay-win;
hi = trueDelay+win;

figure;
maskA = (lags_xcorr>=lo & lags_xcorr<=hi);
plot(lags_xcorr(maskA),r_xcorr(maskA)./max(abs(r_xcorr)),'LineWidth',1); hold on;

maskB = (lag_builtin>=lo & lag_builtin<=hi);
plot(lag_builtin(maskB),r_builtin(maskB)./max(abs(r_builtin)),'LineWidth',1);

maskC = (lags_manual>=lo & lags_manual<=hi);
plot(lags_manual(maskC),r_manual(maskC)./max(abs(r_manual)),'LineWidth',1);

grid on;
xlabel('Lag (samples)');
ylabel('Normalized Correlation');
title('峰值附近放大对比');
legend('xcorr','gccphat','手动GCC-PHAT','Location','best');

结果：