TDOA - GCC-PHAT方法

TDOA（Time Difference of Arrival，到达时间差）是一个在语音前端信号处理中很常见的名词，和它比较相像，也同样重要的还有一个DOA（Direction of Arrival，到达方向）。在麦克风拓扑结构已知的情况下，两者可以相互转换。TDOA表示不同通道之间语音信号的延迟程度，是一个相对的概念，一般用$\tau$表示。最简单的beamforming方法delay and sum中的delay针对的就是它。DOA表示信号源相对麦克风的到达角度，一般用$\theta$表示，如果是线阵的话，范围在0到180度之间。如下图所示，假设麦克风1和麦克风2之间的距离为$\delta$，DOA为$\theta_d$，声速$c$，那么$\tau_{12}$可以有简单的数学计算得到：

$\tau_{12} = \frac{\cos\theta_d \cdot \delta}{c}$

那么DOA/TDOA这个信息怎么用呢？个人觉得可以从下面几个角度理解：

如果目标声源只有一个，外加一些环境噪声，那么TDOA可以用来指导beamformer对该方向的信号进行增强。
如果存在多个说话人，且说话人的方位不同，比如圆桌会议这种，那么：
- 估计的TDOA的数目可以认为是说话人的数量（speaker counting）
- 通过跟踪TDOA的变化，追踪特定的说话人
- 对同一段语音，对估计的DOA方向进行增强，增强说话人

当然现实可能不会这么理想，如果说话人距离较远，说话人方向重合，以及环境复杂（噪声，混响）情景下，都会使得估计方法失效或者beamformer不理想等等。鸡尾酒会问题确实是语音这块一个非常难啃的骨头。现在学术界比较流行的一种Blind的方法就是在不知道声源信息和阵列拓扑的情况下，通过一些先验，或者概率模型估计beamforming中的统计量进行计算，我做过的有CGMM，NN-mask+MVDR/GEV等等，在一些数据集上（比如CHiME4）也取得了很好的效果。

在这里我主要说一个最简单的TDOA估计方法，GCC-PHAT。

首先说一下互相关（CC，Cross-Correlation）。信号处理中，互相关可以用来描述两个信号之间的相似性。离散信号$x_k, y_k$的互相关函数定义为：

$R^{xy}_\tau = E[x_k y_{k + \tau}] = \sum_k x_k y_{k+\tau}$

从上式可知，CC是相对延时$t$的一个函数，显然，可以取使得互相关系数最大的延时值作为TDOA的估计，即：

$\tau = \arg \max_\tau R^{xy}_\tau$

现在普遍认为，这种算法在实际中容易受到噪声和混响的影响，表现不是很稳定，因此引入了广义互相关的概念（GCC，Generalized Cross-Correlation）。和CC一样，TDOA通过获得最大化GCC函数的延迟值得到：

$\mathbf{\tau} = \arg \max_{\tau} \widehat{R}^{xy}_\tau$

其中$\widehat{R}^{xy}_t $定义为generalized cross-spectrum的IDTFT：

$\widehat{R}^{xy}_\tau = \mathcal{F}^{-1}[\gamma_f \cdot \phi^{xy}_f ]$

$\gamma_f$为频域的权重方程，$\phi_f^{xy}$表示cross-spectrum，定义为$E[X_fY_f^*]$，$*$表示共轭。如果将$\gamma_f$定义为$\frac{1}{|\phi_f^{xy}|}$，那么这种估计方法称为为GCC-PHAT。PAHT表示phase transform，因为做了幅值的归一化之后，相当于只留下了相位信息。

语音信号处理中，我们通常在短时频域上进行。在每一帧上，进行一次$\widehat{R}^{xy}_\tau$计算，可以得到所谓的angular spectrogram：

$\varPsi_{\tau t} = \sum_f \overline{H_{tf}} e^{j2\pi f\tau}$

其中$H_{tf} = X_{tf} \odot Y_{tf}^*$，$\overline{H_{tf}}$表示做幅值归一化。实际计算中，构建矩阵$\varOmega \in R^{\varTheta \times F}$，使得

$\varOmega_{\tau f} = e^{j2\pi f\tau}$

之后计算矩阵乘法$\varPsi =\overline{H} \varOmega$即可。下面贴一段我写的可视化angular spectrogram的程序，旨在说明计算逻辑：

#!/usr/bin/env python
# coding=utf-8
# wujian@2018

import argparse
import glob

import numpy as np
import matplotlib.pyplot as plt
import utils

speech_speed    = 340.29
distance_of_mic = 1
samp_frequency  = 16000

def compute_stft(wave_path, frame_length, frame_shift):
    wave_samp = utils.audio_read(wave_path)
    return utils.stft(wave_samp, frame_length, frame_shift)

def compute_angspect(coherence, num_tdoa):
    """
        coherence: num_bins x num_frames
        return num_tdoa x num_frames
    """
    # num_bins x num_toda
    num_bins, num_frames  = coherence.shape
    max_tdoa = distance_of_mic / speech_speed
    # from 0 to 180
    est_tdoa = np.linspace(-max_tdoa, max_tdoa, num_tdoa)
    sep_freq = np.linspace(0, samp_frequency / 2, num_bins)
    # num_tdoa x num_bins
    exp_part = np.outer(2j * np.pi * est_tdoa, sep_freq)
    return np.dot(np.exp(exp_part), coherence).real

def run(args):
    stft = [compute_stft(wave, args.frame_length, args.frame_shift) for wave in sorted(glob.glob(args.pattern))]
    assert len(stft) == 2, 'Check pattern {} please'.format(args.pattern)
    coherence = stft[0] * stft[1].conj() / (np.abs(stft[0]) * np.abs(stft[1]))
    angspect = compute_angspect(coherence， args.num_tdoa)
    angspect[angspect < 0] = 0
    plt.subplot(211)
    plt.ylabel('TDOA index')
    plt.xlabel('Frame index')
    plt.title('GCC-PATH Angular Spectrogram')
    plt.imshow(angspect, cmap=plt.cm.binary, aspect='auto', origin='lower')
    plt.subplot(212)
    plt.ylabel('Frequency index')
    plt.xlabel('Frame index')
    plt.title('Log-Magnitude Spectrogram')
    plt.imshow(np.log(np.abs(stft[0])), cmap=plt.cm.jet, aspect='auto', origin='lower')
    plt.show()


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description="Command to compute angular spectrogram using GCC-PHAT methods")
    parser.add_argument('pattern', type=str, help='location of 2-channel wave files')
    parser.add_argument('--num-of-tdoa', type=int, default=128, dest="num_toda",
                        help='Number of todas to estimate')
    parser.add_argument('--frame-shift', type=int, default=256, dest="frame_shift",
                        help='Frame shift in number of samples')
    parser.add_argument('--frame_length', type=int, default=1024, dest="frame_length",
                        help='Frame length in number of samples')
    args = parser.parse_args()
    run(args)

注：程序无法完整运行，utils.py我未给出。

在一些比较理想或者simulate的数据下，这种方法work的很好，这部分例子可以参见“Blind Speech Separation and Enhancement With GCC-NMF ”这系列论文以及github。实际的数据中，表现还是差强人意的，我在实录的数据上只在部分条件下（比如说话人距离麦克风较近等等）可以提供较为准确的TDOA估计信息。比如下面几种情况：

上图中的angular spectrum中在index-80的位置，在时间轴上出现几段较为明显的峰值，由此可以得到以下信息：

峰值出现的所在的index即为估计的最佳TDOA索引。
峰值不连续处表示该说话人尚未说话（起到类似VAD的作用）。
峰值所对应的index分布向着90处移动，说明该说话人可能相对麦克风在移动。

如果在不同index处，同时出现峰值，那么就说明存在多个说话人，比如下面这种情况，可以认为在这段语音中，可能存在三个说话人。

GCC-PHAT方法的局限性还表现在仅仅只能依据两路信号进行估计。如果现在有多路麦克风，那么就需要在原有算法的基础上进行改进。这部分由于我还没有具体的实践，暂时无法继续深入。

之前也提到过，Blind的方法目前已经摆脱对具体TDOA/DOA的依赖的，以MVDR beamformer为例子，往往会用noise的协相关矩阵的最大特征值作为steer vector的估计。而在传统的信号处理中，TDOA是形成steer vector不可或缺的一个因素，假如已经计算出相对延迟序列$\tau_1, \tau_2, \cdots, \tau_{N-1}$，那么steer vector为：

$\mathbf{d}_{\theta} = [1, e^{j\omega \tau_1}, e^{j\omega \tau_2}, \cdots, e^{j\omega \tau_{N-1}}]$

在beamformer设计中很重要的一个东西是beam pattern，在我们做Blind Source Separation/Enhancement的时候，往往会忽视这个分析手段，只看重最终的WER，在后面我会开一篇文章，介绍使用TDOA进行beamformer，以及使用beam pattern进行分析的方法。什么时候更呢……等哪天准备工作做得充分了吧。