---

前言

在实际工程中，很多人在实现 麦克风阵列声源定位（DOA） 时都会遇到一个问题：

角度输出不稳定，甚至出现明显跳动。例如在声源基本不移动的情况下，系统输出的角度却可能在：

45° → 70° → 50° → 65°

之间不断变化。

一开始很多人会认为是 GCC 算法实现问题，但实际上问题往往来自多个方面：

例如：

• 麦克风阵列几何结构
• TDOA 延迟量化误差
• 语音信号非平稳性
• 噪声和混响

在上一篇文章四麦克风声源定位实战：基于 GCC-PHAT + 最小二乘法实现 DOA中，详细介绍了基于四麦实现的TDOA的原理和代码。

在本文中，我们在相同的实验条件下，对 三麦克风阵列和四麦克风阵列 的定位效果进行了对比实验，并分析导致 DOA 跳动的几个关键原因：

• 阵列几何约束
• TDOA 量化误差
• 算法与环境影响

通过实验与理论分析，可以更直观地理解：为什么三麦克风定位更容易出现角度跳动，以及工程中应该如何优化。
|版本声明：山河君，未经博主允许，禁止转载

---

一、三麦与四麦 DOA 实验对比

1.实验

为了分析 DOA 跳动问题，我们在相同实验条件下，对 三麦克风阵列和四麦克风阵列 的定位结果进行了对比。

实验条件如下：

• 麦克阵列：四麦为{ 0.00f, 0.00f},{ 0.02f, 0.00f},{ 0.00f, 0.02f}, { 0.02f, 0.02f}，三麦去除{ 0.02f, 0.02f}
• 采样率：48k
• 输入信号：同一段语音，三麦去除对应声道
• 每1024帧输出一个角度

2.三麦定位公式

$$\begin{gathered} d_{12}=(x_2-x_1)\cos \theta +(y_2-y_1)\sin \theta \\ {d}_{13}=(x_3-x_1)\cos \theta +(y_3-y_1)\sin \theta \end{gathered}$$
写成矩阵：
$$\begin{gathered} \left[\begin{array}{c}{x}_2-x_1,y_2-y1\\ x_3-x_1,y_3-y1\end{array} \\ \right]\left[\begin{array}{c}\cos \theta \\ \sin \theta \end{array} \\ \right]=\left[\begin{array}{c}{d}_{12}\\ d_{13}\end{array}\right] \end{gathered}$$

设
$$A=\left[\begin{array}{c}{x}_2-x_1,y_2-y1\\ x_3-x_1,y_3-y1\end{array}\right]d=\left[\begin{array}{c}{d}_{12}\\ d_{13}\end{array}\right]$$
那么
$$\begin{gathered} \left[\begin{array}{c}\cos \theta \\ \sin \theta \end{array} \\ \right]=A^{-1}d \end{gathered}$$

3.对应代码

1）矩阵求逆

// 2x2矩阵求逆
    float dx1 = m_mic[1].x - m_mic[0].x;
    float dy1 = m_mic[1].y - m_mic[0].y;
    float dx2 = m_mic[2].x - m_mic[0].x;
    float dy2 = m_mic[2].y - m_mic[0].y;
    float det = dx1 * dy2 - dx2 * dy1;
    m_finv11 =  dy2 / det;
    m_finv12 = -dy1 / det;
    m_finv21 = -dx2 / det;
    m_finv22 =  dx1 / det;

2）计算角度

// A^{-1}d
float DoaEstimator3Mic::solveAngle(const std::vector<float>& tau)
{
    float b1 = SOUND_SPEED * tau[0];
    float b2 = SOUND_SPEED * tau[1];


    float cos_theta = m_finv11 * b1 + m_finv12 * b2;
    float sin_theta = m_finv21 * b1 + m_finv22 * b2;

    // 归一化（非常重要）
    float norm = sqrtf(cos_theta * cos_theta + sin_theta * sin_theta);
    if (norm > 1e-6f)
    {
        cos_theta /= norm;
        sin_theta /= norm;
    }

    return atan2f(sin_theta, cos_theta);
}

4.输出结果比较

上图中左图是三麦阵列，右图是四麦阵列，输入是同一段语音，由输出结果对比可以发现：

• 两种阵列都存在角度跳动
• 三麦阵列的角度波动明显大于四麦阵列

而出现这些问题的原因总的来说有几下几点原因：

• 几何约束：受到物理影响
  • 阵列数量
  • 阵列位置
• TDOA 量化：受到策略影响
  • 采样率
  • 帧率
• 算法选择：受到不同GCC算法权重影响
  • 噪声敏感
  • 混响敏感

二、几何约束

1.阵列数量选择

对于四麦DOA为超定方程：
$$\begin{gathered} A_{3\cdot 2}\left[\begin{array}{c}\cos \theta \\ \sin \theta \end{array} \\ \right]=d=>u=(A^{T}A{)}^{-1}{A}^{T}b \end{gathered}$$
其优势：

• 具有冗余约束
• 可以抑制 TDOA 噪声

而三麦方程数等于未知数，只存在唯一解，所以任意TDOA误差都会直接影响结果，造成角度跳动明显

2.位置排布

1）阵列分布

麦克风的排布应该遵循：

• 基线方向越均匀越好
• 阵列尽量对称

例如：四麦克往往是正方形，而三麦最优是等边三角形，等边三角形阵列具有最好的方向均匀性。

上文示例中的三麦排列为直角三角形，就导致两个基线接近正交，某些方向几何敏感度差。

2）阵列间距

麦克直接的间距不能过大或者过小。例如：

• 过小：间距1cm，最大延迟$\tau =d/c$，对应48k采样有效采样点为1~2sample
• 过大：间距10cm，$f>c/2d$，会出现空间混叠，语音高频会产生多解。

工程经验（语音阵列）：2 cm ~ 5 cm

三、TDOA 量化

1.精度误差

1）误差分析

假设采样率为48k，那么：

• 每个采样点代表的延迟$T_s=1/48000=20.8\mu s$
• 声速$c=343m/s$
• 一个采样点对应的传播距离$\Delta d=c{T}_s=343\cdot 20.8\cdot 10^{-6}\approx 7.1mm$

假设系统分辨的延迟以及对应的距离差为

... 
-2 sample -> -14mm
-1 sample -> -7mm
0 sample -> 0mm
1 sample -> 7mm
2 sample -> 14mm
...

如果真实误差是在10mm，则会被量化到最近的 sample，也就是最接近的7mm或者14mm，根据计算角度公式：
$$\theta =\arccos (d/D)$$
算出的角度为：

• 真实：$\theta =\arccos (10/20)=60^{\circ }$
• 7mm：$\theta =\arccos (7/20)=69.5^{\circ }$
• 14mm：$\theta =\arccos (7/20)=45.6^{\circ }$

即看到的输出角度很有可能在$69.5^{\circ }\to 45.6^{\circ }$跳动。

2）解决策略：亚采样插值

抛物线插值的核心原理：假设相关函数峰值附近可以用一个二次函数近似，即
$$R(x)\approx a{x}^2+bx+c$$
如果我们已知：

• 峰值点：$R(0)$
• 左边一点：$R(-1)$
• 右边一点：$R(+1)$

就可以拟合出一个抛物线，然后求它的顶点位置。例如：

y1 = R(k-1)
y2 = R(k)
y3 = R(k+1)

那么抛物线顶点偏移量：
$$\delta =\frac{y_1-y_3}{2(y_1-2y_2+y_3)}$$
最终亚采样后的位置在：
$$\tau =k+\delta$$
对应代码：

{
....
	int k = max_idx;
    int k_minus = (k - 1 + FRAME_SAMPLE) % FRAME_SAMPLE;
    int k_plus  = (k + 1) % FRAME_SAMPLE;

    float delta = parabolicInterp(m_pDataR[k_minus], m_pDataR[k], m_pDataR[k_plus]);
    float peak = k + delta;
....
}
float DoaEstimator3Mic::parabolicInterp(float y1, float y2, float y3)
{
    float denom = 2.0f * (y1 - 2.0f * y2 + y3);

    if (fabs(denom) < 1e-6f)
        return 0.0f;

    float delta = (y1 - y3) / denom;

    return delta;
}

2.帧率过高

1)误差分析

如果每 1024 点输出一次 DOA，48kHz 下：$1024/48000=21ms$，意味着每秒输出约 47 个角度。而语音信号具有：

• 非平稳性
• 噪声
• 混响

导致每帧的 TDOA 会波动。 从而出现： 角度快速跳动。

为了对抗这种现象，需要使用者首先理解自己是处于什么样的场景下。很多工程系统不会每帧输出 DOA， 而是 200ms 或 500ms 输出一次稳定结果。

2）直方图策略

直方图策略的核心是针对于帧数统计，如果使用是在一个稍长时间快速定位，例如快速$A⟶B$

思路：

• 建立 0~360° 的 bin（比如 5° 一格）
• 统计最近 N 帧 DOA
• 找出现次数最多的 bin
• 输出该 bin 的中心角

这相当于一个：

• 模式滤波器（mode filter）
• 比均值滤波强，因为均值会被离群点拉偏。

代码示例：

#include "angle_histogram.h"
#include <fftw3.h>
#include <math.h>

AngleHistogram::AngleHistogram(int nWindow, float fBinWidth)
    :m_window(nWindow),
    m_binWidth(fBinWidth),
    m_binCount(static_cast<int>(360.0f / fBinWidth))
{
    m_hist.resize(m_binCount);
    std::fill(m_hist.begin(), m_hist.end(), 0);
}

float AngleHistogram::update(float fAngle)
{
    int newIdx = angleToBin(fAngle);
    m_hist[newIdx]++;
    m_buffer.push_back(newIdx);

    if (m_buffer.size() > m_window)
    {
        int oldIdx = m_buffer.front();
        m_buffer.erase(m_buffer.begin());
        m_hist[oldIdx]--;   // 减去旧bin
    }

    int maxIdx = findMaxBin();
    return binCenter(maxIdx);
}

int AngleHistogram::findMaxBin()
{
    int max_index = 0;
    if (!m_hist.empty()) {
        auto max_it = std::max_element(m_hist.begin(), m_hist.end());
        max_index = std::distance(m_hist.begin(), max_it);
    }

    return max_index;
}

float AngleHistogram::binCenter(int index)
{
    return (index + 0.5f) * m_binWidth;
}

结果：

3）指数平均

对于想短时定位，或者想看到$A⟶B$的轨迹，则可以进行指数平均，即对于上一帧的角度和当前帧角度施加不同的权重

float theta = solveAngle(tau);
    float cos_new = cosf(theta);
    float sin_new = sinf(theta);
    m_fCos = (1 - ALPHA) * m_fCos + ALPHA * cos_new;
    m_fSin = (1 - ALPHA) * m_fSin + ALPHA * sin_new;
    return atan2f(m_fSin, m_fCos) * 180.0f / M_PI;

结果：

4）圆形平均

这种方式就介于指数平均和直方图之间了，也就是介于平滑轨迹和快速定位之间。其原理是

把角度看成单位圆上的点：

x = cos(theta)
y = sin(theta)

累计：

X = Σ cos(theta_i)
Y = Σ sin(theta_i)

最后角度：

theta = atan2(Y, X)

其角度变化大的相当于权重更大。

其实现示例：

#ifndef CIRCULARMEAN_H
#define CIRCULARMEAN_H

#include <deque>
#include <cmath>

class CircularMean
{
public:

    CircularMean(int window)
        :m_window(window), m_sumX(0), m_sumY(0)
    {}

    float update(float angleDeg)
    {
        float rad = angleDeg * M_PI / 180.0f;

        float x = cos(rad);
        float y = sin(rad);

        m_buffer.push_back({x,y});

        m_sumX += x;
        m_sumY += y;

        if (m_buffer.size() > m_window)
        {
            auto old = m_buffer.front();
            m_buffer.pop_front();

            m_sumX -= old.first;
            m_sumY -= old.second;
        }

        float angle = atan2(m_sumY, m_sumX) * 180.0f / M_PI;

        if (angle < 0)
            angle += 360;

        return angle;
    }

private:

    int m_window;

    float m_sumX;
    float m_sumY;

    std::deque<std::pair<float,float>> m_buffer;
};

#endif // CIRCULARMEAN_H

结果：

5）VAD检测

以上所有的策略都建议增加：结合 VAD 门控

if (VAD == active)
    push angle to histogram
else
    不更新

否则静音帧会污染角度，这里可以参考webrtc vad实现方式：webrtc之语音活动下——VAD人声判定原理以及源码详解

四、算法误差

1.误差分析

1）麦克一致性误差

真实麦克风存在：

• 增益误差
• 相位误差
• 频响差异

这些误差会影响：互相关函数形状，从而影响 TDOA。

2）多路径/非平稳（混响/噪声）

室内环境通常存在：

• 墙面反射
• 桌面反射
• 天花板反射

或者

• 间歇信号
• 有停顿
• 有爆破音
• 低SNR

因此 GCC 会出现，多个相关峰，例如：

• 主峰 (直达声)
• 副峰1、2、3、4
• 最大峰会跳到副峰
• 互相关峰值会变得模糊
• 多个峰值

2.解决方案

选择合适的GCC算法：声源定位核心算法：GCC-PHAT / SCOT / Roth / ML 原理与对比。

---

总结：

通过三麦克风与四麦克风阵列的实验对比，可以发现 DOA 跳动问题通常来自多个方面的叠加影响，主要包括：

• 阵列几何约束
• TDOA 量化误差
• 语音信号非平稳
• 噪声与混响

工程建议，在语音阵列定位系统设计中，通常建议：

• 麦克风数量 ≥ 4
• 麦克间距 2 cm ~ 5 cm
• 使用 GCC-PHAT + 亚采样插值
• 增加 时间稳定策略
• 配合 VAD 门控

这样可以显著提升 DOA 输出的稳定性。

如果本文对你有所帮助，欢迎点赞、收藏和关注。
后续会继续分享麦克风阵列声源定位相关内容。