初识声源识别-TDOA及GCC-PHAT简述

本文主要介绍了一些声源识别的基础概念，常用方法，以及传统的运用广义互相关GCC的TDOA方法。

一、相关定义

混响

• 定义：声波在室内传播时，被障碍物反射，每一次被反射都要被障碍物吸收一些。当声源停止发声后，声波要在室内经过多次反射和吸收后才会消失。并且在声源停止发声后，还有若干个声波混合持续一段时间（即声延续现象）。这种现象叫做混响。

相关函数

• 定义：用来描述信号X(t)，Y(t)（这两个信号可以是随机的，也可以是确定的）在任意不同两个时刻s、t的取值之间的相关程度。两个信号之间的相似性大小用相关系数来衡量：

  • \(\rho = \frac{COV(X,Y)}{\sqrt D(X)\sqrt D(Y)}=\frac{\sigma_{xy}}{\sigma_x \sigma_y}\)
  • 若相关系数=0，则代表X与Y不相关，相关性系数越大，相关性也就越大。

• 自相关函数：描述随机信号x(t)在任意不同时刻t₁，t₂的取值之间的相关程度。是信号自身的互相关，表示同一序列在不同时刻的相关程度，是用来寻找重复模式的数字工具。常用作信号处理中的分析函数。

  • 定义式：\(R(s,t) = E(X(s)*X(t))\)
  • 或者:\(R_{xx}(\tau)=\lim\limits_{T\rightarrow \infty}\frac{1}{T}\int^T_0x(t)x(t+\tau)dt\)
    • 对于周期信号，积分平均时间T为信号周期。在有限时间内，单个脉冲在T趋于无穷大的情况下，该平均值将无限接近于0。此时自相关函数可以表达为\(R_xx(\tau)=\int^{+\infty}_{-\infty} x(t)x(t+\tau)dt\)。
    • 自相关函数就是信号x(t)和它的时移信号\(x(t+\tau)\)的乘积的平均值

• 互相关函数：描述随机信号x(t)，y(t)在任意两个不同的时刻s，t取值之间的相关程度。

  • 定义式：\(R(s,t)=E(X(s)*Y(t))\)
  • 连续函数：\((f * g)(\tau)=\int^{+\infty}_{-\infty}f^*(t)g(t+\tau)dt\)
  • 离散函数：\((f*g)(n)=\sum\limits_{-\infty}^{+\infty}f^*[m]g(n+m)\)

近场模型

• 将声波看作球面波，主要考虑麦克风阵列各阵元接收到信号的幅度差。
• 主要考虑位置关系

远场模型

• 将声波看作平面波，忽略阵元接收信号间的幅差，近似认为各阵元接收信号之间为简单的时延关系。
• 主要考虑时间关系

时域

• 时间域，自变量是时间，即横轴是时间，纵轴是信号的变化。其动态信号是描述信号在不同时刻取值的函数。时域分析是以时间轴为坐标表示动态信号关系。

频域

• 频率域，自变量是频率，即横轴是频率，纵轴是该频率信号的幅度，该图就是常说的频谱图。频谱图描述了信号的频率结构以及频率与该频率信号幅度的关系。频域是把时域波形的表达式作傅立叶变换得到复频域的表达式，所画出的图就是频谱图。

相位谱

• 相位是信号组成的重要部分。相位在信号图中表达为时域图像中正弦波离频率轴最近的波峰在底面的投影到轴的距离。
• 将投影得到的时间差依次除以除以所在频率的周期，可以得到相位谱。
• 相位谱中的相位除了0就是Pi，因为cos(t+Pi)=-cos(t)。因此相位为Pi的波是上下翻转了的。

傅立叶级数

• 傅立叶级数是一种对周期信号进行分解的方式。通过对信号进行分解，可以更直观和方便的对某段信号进行分析。
• 傅立叶级数在时域是一个周期的且连续的函数，在频域是一个非周期的且离散的函数
• 傅立叶把任何一个波都看成是数个正弦波的叠加

欧拉公式

• \(e^{it}=cos(t)+i*sin(t)\)将正弦波统一成了简单的指数形式，它描绘的是一个随着时间变化，在复平面上做圆周运动的点，在时间轴上形成了一条逆时针旋转的螺旋线。其实部是一个余弦函数，虚部部分是一个正弦函数。
• 借助欧拉公式及其另一种形式\(e^{-it}=cos(t)-i*sin(t)\)，我们可以获得\(\cos(t) = \frac{e^{it}+e{-it}}{2}\)。
• 逆时针旋转的称为正频率，顺时针旋转的称为负频率

傅立叶变换

• 傅立叶变是频率从0到无穷的叠加，是一种从时间转换为频率的变化。
• 连续傅立叶变换
  • 表示将平方可积的函数f(t)表示成复指数函数的积分或级数形式
  • \(F(\omega)=\int\limits^{\infty}_{-\infty}f(t)e^{-i\omega t}d\omega\)将频率域的函数\(F(\omega)\)表示成时间域的函数f(t)的积分形式
  • 连续傅立叶变换的逆变换为\(f(t)=\frac{1}{2\pi}\int\limits^{\infty}_{-\infty}F(\omega)e^{i\omega t}d\omega\)。将时间域函数f(t)表示为函数\(F(\omega)\)的积分
  • 连续傅立叶变换是傅立叶级数的推广
• 离散时域傅立叶变换
  • DTFT在时域上离散，在频域上是周期的。DTFT可以被看作是傅立叶级数的逆变换
• 离散傅立叶变换
  • DFT是连续傅立叶变换在时域和频域上都是离散的形式，将时域信号上的采样变换为在离散时间傅立叶变换频域上的采样
• 傅立叶变换分类
  • 根据信号的不同类型，可以把傅立叶变换分为四种类型
    1. 非周期性连续信号 傅立叶变换
    2. 周期性连续信号 傅立叶级数
    3. 非周期性离散信号 离散时域傅立叶变换
    4. 周期性离散信号 离散傅立叶变换

互相关

• 对两个函数分别作复数共轭和反向平移并使其相乘的无穷积分，或者说：第一个函数依次作复共轭和平移后与第二个函数相乘的无穷积分。用以比较两个信号的相似程度

• 在更通俗的说法中，互相关就像是在两段信号之间滑动一个窗口，然后看在不同时间点窗口下两个信号的形状有多相似。如果窗口下的信号形状非常相似，互相关值就会很高，表示两个信号在那个时间点上很相像。

  在音频处理中。例如，如果有两个麦克风录制的声音，就可以使用互相关来找出声音中的时间延迟。

  总的来说，互相关是一种比较两个信号在不同时间点上相似程度的方法。

• 计算互相关函数的步骤

  1. 信号准备：将两个信号进行准备，确保在时间上对齐。
  2. 傅立叶变换：对准备好的信号进行傅立叶变换，将他们转换到频域。这一步通常使用快速傅立叶变换算法(FFT)实现。
  3. 复共轭：对一个信号进行复共轭操作，保持实部不变，虚部取反
  4. 点乘：将第一个信号与第二个信号的复共轭进行点乘。等价于在频域上将幅度相乘，相位相加。
  5. 逆傅立叶变换：将点乘后获得的频域上的信号转换为时域。
  6. 归一化：对第5步中的结果进行归一化处理，确保互相关函数的值在合适的范围内。

• 互相关系数：代表了两个信号之间的相似性或相关性，用以衡量两个信号在时间上的相对偏移。

  • 互相关系数的计算涉及对两个信号进行点乘或卷积运算。它的值范围在 -1 到 1 之间，其中：
    • 当互相关系数接近 1 时，表示两个信号具有高度正相关性，即它们在时间上趋向于同时增大或减小。
    • 当互相关系数接近 -1 时，表示两个信号具有高度负相关性，即它们在时间上趋向于一个增大时，另一个减小。
    • 当互相关系数接近 0 时，表示两个信号之间没有线性相关性。
  • 表达式
    • \(\rho_{xy} = \frac{\sum_{t=1}^{N}(x_t - \bar{x})(y_{t+\tau} - \bar{y})}{\sqrt{\sum_{t=1}^{N}(x_t - \bar{x})^2 \sum_{t=1}^{N}(y_{t+\tau} - \bar{y})^2}}\)
    • (\(\rho_{xy}\)) 表示互相关系数，(\(x_t\)) 和 (\(y_t\)) 分别表示两个信号在时间点 (t) 的取值，(\(\bar{x}\)) 和 (\(\bar{y}\)) 分别表示两个信号的均值，(\(\tau\)) 表示时间偏移，(N) 是信号的长度。这个公式表示的是标准的互相关系数的计算方式。

广义互相关

• 在互相关的基础上位了减少噪声和混响带来的影响，在互功率谱使用频域加权函数进行计算，然后通过IFFT运算找到峰值，峰值即时延
• 
• 公式
  • GCC数学表达式：\(R_{x_1x_2}(\tau)=\int^{\infty}_{-\infty}\psi_{12}(\omega)G_{x_1x_2}(\omega)e^{j\omega \tau}d\omega\)

CCF

• 交叉相关函数，用于衡量两个信号之间在不同空间或时间点上的相似性或相关性
• 如果 (\(CCF(\tau)\)) 在某个时间延迟 (\(\tau\)) 处有较大的正值，表示两个信号在该时间延迟下存在正相关性；如果有较大的负值，表示存在负相关性。零值表示没有明显的相关性。

二、声源定位

1.声源定位的定义

• 根据由按照一定几何拓扑结构放置的麦克风阵列记录的多通道声学信号估计一个或者多个声源相对于某个任意参考位置的位置问题，即对声源位置的估计。
• 通常情况下被简化为声源到达方向(DOA)的估计，即专注于方位角和仰角，不估计声源到达麦克风阵列的距离。

2.声源定位的传统方法

声源定位的传统方法都是基于信号/信道模型和信号处理技术，但是传统方法在噪声，混响和多个发声源的情况下表现不佳。

TDOA：

• 已知麦克风阵列的几何形状，根据时延估计(声学信号到达麦克风阵列中不同麦克风之间的时间差)来进行DOA估计。
• 步骤：
  1. 时延估计
     • 通过将阵列中麦克风两两组合，使用互相关，广义互相关和相位差等方法计算出传播延迟。
     • 常用时延估计法：
       • 广义互相关(GCC)
       • 互功率谱相位法(CSP)
       • 特征值分解
       • 广义特征值分解
       • 声学传递函数比
  2. 定位估计
     • 根据麦克风阵列的几何结构和第一步已知的时延信息，通过搜索或者几何方法，来进行声源的方向和距离估计。
  3. 特点：受到噪声，混响，多径等因素影响，导致定位不准确
  4. GCC-PHAT：
     • 具有相变的广义互相关方法。若用于2-麦克风序列，则该值为两个麦克风接收信号之间的交叉功率谱(CPS)的加权逆傅立叶变换。通过找到使GCC-PHAT函数值最大的时延来获得TDOA估计。

BF：

• 基于可控波束形成的声源定位方法，通过对麦克风阵列采集到的数据进行滤波，加权，叠加，获得一个具有指向性的波束。
• 可以通过方法某个特定方向上的数据来减少其他方向上的干扰。
• DSB算法：
  • 通过对不同麦克风信号之间的相对延迟进行补偿，将叠加延时后的信号形成一个输出。
  • 特点：该方法可以有效减少非相关噪声影响，增强目标方向信号强度，提高信噪比。
  • 还有更复杂的可控波束形成起，可以同时对信号进行延时补偿和滤波处理。
• 特点：采用遍历式搜索，计算代价较大，实际应用中难以保证实时性。

MUSIC：

• 基于高分辨率谱估计的声源定位方法，分析每个元素信号之间的相关性，拟合出一个矩阵，通过矩阵来确定声源位置。
• 特点：克服了空间分辨率受限问题，但是在极短时间内采集满足要求的数据较为困难。

SRP：

• 构建声学信号功率图，通过指向延迟和总和波束形成器，朝向标准网格中的候选位置，测量在这些位置上的能量(声音强度)，其中的极大值就是主要对应的源的DOA估计。
• 特点：计算代价高，类似于网格搜索。

子空间方法：

• 是一种定位算法类，类中的方法依赖于基于平均时间的CPS矩阵。

  • MUSIC(多信号分类)方法:

  • 应用EVD(特征值分解)来估计信号和噪声子空间，然后使用信号和噪声子空间来估计指定方向上是否存在源，即使用空间滤波波束形成。

  • ESPRIT算法

  • 通过松弛信号参数，利用源的子空间结构来直接推断源DOA。

• 特点：对噪声的具有鲁棒性(健壮性)，但是对混响敏感

GMM(基于概率生成混合模型)方法：

• 该模型是一个高斯混合模型的变体，每个源都有一个高斯分量来进行定位，筛选源的位置。
• 期望最大化(EM)算法：
  • 利用时间频率域中声源的稀疏性。
• GMR算法：
  • 一种直接以回归模式运行的GMM变体。
  • 是局部线性，全局非线性。
• 该方法与贝叶斯推理息息相关。

ICA：

• 通过假设和利用他们的相互统计独立性来检索包含混合的异源信号。
• 常用于对音频处理以进行盲源分离，也对多源SSL问题有用。

三、TDOA及GCC-PHAT

1.TDOA的声源定位方法及原理

假设在空间中有声源记为s(t)，其在空间位置为S。两个麦克风m₁,m₂，在空间中位置分别为M₁,M₂，接收到的信号为x₁(t),x₂(t)

2.GCC-PHAT算法数学模型

远场模型示意图如下:

声音信号s(k)位于远场，视为平面波，入射角为\(\theta\)，两个麦克风y₁(k),y₂(k)之间的距离为d

时延表达式:\(\tau_{12}=d\cos{\theta}/c\)

计算出入射角\(\theta\)就可以定位出声音信号的位置。

近场声源可以估计声源到达传感器的角度，并且可以估计声源与麦克风之间的距离

远场声源可以估计入射角度，但是难以确定声源与麦克风阵列之间的距离

1)互相关方法

先定义两个观测信号y₁(k)与y₂(k)之间的互相关函数

• \(R^{cc}_{y_1y_2}(p)=E[y_1(k)y_2(k+p)]\)
• 当p=\(\tau\)时，\(\tau\)为相对时延，互相关值达到最大
  • \(\hat{\tau}^{CC}=\arg\mathop{\min}\limits_{p}r^{CC}_{y_1y_2}(p)\)
• 互相关方法易受到噪声，混响等的影响

2)广义互相关方法

再定义广义互相关方法，此时两个麦克风之间的TDOA估计可以等效为能够使麦克风输出的滤波信号之间的CCF最大的时间间隔：

• \(\hat{\tau}^{GCC}=\arg\mathop{\max}\limits_{\tau}r^{CC}_{y_1y_2}(p)\)

GCC函数:

• \(r^{GCC}_{y_1y_2}(p)=F^{-1}[\psi_{y_1y_2}(f)]=\int^{\infty}_{-\infty}\psi_{y_1y_2}(f)e^{j2\pi fp}df=\int^{\infty}_{-\infty}\vartheta(f)\phi_{y_1y_2}(f)e^{j2\pi fp}df\)
• 广义互频谱：\(\psi_{y_1y_2}(f)=\vartheta(f)\phi_{y_1y_2}(f)\)
• 互频谱：\(\phi_{y_1y_2}(f)=E[Y_1(f)Y_2^*(f)]\)
• 这个公式是cc函数的互相关函数傅立叶变换

3）相位变换

• TDOA的估计是通过互频谱相位而不是幅度，所以要简单舍弃幅度，保留相位
• \(\vartheta(f)=\frac{1}{|\phi_{y_1y_2}(f)|}\)
• 理想GCC函数:
  • \(r^{\text{phat}}*{y_1y_2}(p) = \int^{\infty}_{-\infty} e^{-j2\pi f(p-\tau)} \delta(p - \tau) , df = \begin{cases}\infty, p=\tau\\0, \text{其他}\end{cases}\)