语音数字信号处理——GCC-PHAT计算回声时延

---

一、回声估计的重要性

在回声消除真正处理时，首先就要进行时延估计，这一方便做的比较好的是webrtc，然而类似于speex、硬件3a上并没有对于时延进行估计，这也是为什么实际使用两个相同的pcm文件进行回声消除效果非常的好，但是一但将其中一个pcm文件进行延迟个10ms以上，再进行回声消除，基本上是不起到效果的原因。所以在音视频通话场景下，就需要用到回声估计。

在实际播放的音频信号经过扬声器、空间传播、麦克风拾取，会有一定的时延，这个时候如果直接进行回声消除，可能无法正确匹配回声信号，导致回声建模不正确，或者消除掉非回声部分。所以这个时候就需要对于时延进行估计。

本章中使用GCC-PHAT算法进行回声延迟估计，这种算法只做理论研究，在实际场景中，并不能简单直接使用该算法，因为还要考虑到噪声干扰和回声强度变化导致无法及时进行估计。

|版本声明：山河君，未经博主允许，禁止转载

---

一、设计思路

1.整体思路

计算两个信号的时延，GCC-PHAT利用的是广义互相关-相位变换，计算两个信号的互功率谱，从而得到它们在不同频率分量上的相互关系。

解释一下什么是互功率谱：两个信号在频域中的相关性度量，表示它们在不同频率分量上的相互关系。表达公式为：$R_{xy}=X(f)Y^{\ast }(f)$。
在音频进阶学习九——离散时间傅里叶变换DTFT中知道，对于两个复数如果正交，那么它们彼此垂直或互不相关，此时我们得到了频谱的相关性，但是我们关注的是此时的相位信息，所以对于取共轭的原因是为了计算相位差（不取共轭是相位和）：
$$R_{xy}=X(f)Y^{\ast }(f)=|X(f)||Y(f)|(j{e}^{{\omega }_x-{\omega }_y})$$
此时我们就保留了频谱相同时延迟的信息，这个时候我们去除幅频响应，只保留相位信息，那么还原成时域图像后，时域最大时即为延迟最大时。

其中实现过程中需要使用到傅里叶变换库，这里不做编译和FFT的API解释了。

1.近端信号和参考信号

输入浮点型信号的目的是为了进行三角运算后精度丢失，所以对于整形pcm数据，需要转为浮点型进行运算

for (size_t i = 0; i < vecRef.size(); ++i) {
            doubleRef[i] = static_cast<double>(vecRef[i]);
        }
for (size_t i = 0; i < vecNear.size(); ++i) {
            doubleNear[i] = static_cast<double>(vecNear[i]);
        }

2.进行傅里叶变换

将参考信号和近端信号进行快速傅里叶变换

    const int N = x.size();
    const int N_fft = N / 2 + 1;
    fftw_complex *X = fftw_alloc_complex(N_fft);
    fftw_complex *Y = fftw_alloc_complex(N_fft);
    fftw_plan planX = fftw_plan_dft_r2c_1d(N, const_cast<double*>(x.data()), X, FFTW_ESTIMATE);
    fftw_plan planY = fftw_plan_dft_r2c_1d(N, const_cast<double*>(y.data()), Y, FFTW_ESTIMATE);
    fftw_execute(planX);
    fftw_execute(planY);

其中const int N_fft = N / 2 + 1，是由于FFT计算时，由于共轭对称性，FFT的计算结果是对称的，所以FFT的计算结果只需要保存前半部分，而后半部分作为冗余信息进行节省。

至于FFT的原理与设计，可以查看音频进阶学习二十一——FFT（DIT-FFT以及实现代码）这篇文章。

3.计算互功率谱

计算参考信号和近端信号的互功率谱

    fftw_complex *R = fftw_alloc_complex(N_fft);
    for (int i = 0; i < N_fft; ++i) {
        // 计算互功率谱
        const double X_real = X[i][0];
        const double X_imag = X[i][1];
        const double Y_real = Y[i][0];
        const double Y_imag = Y[i][1];

        // X * conj(Y)
        const double cross_real = X_real * Y_real + X_imag * Y_imag;
        const double cross_imag = X_imag * Y_real - X_real * Y_imag;
    }

其中X_real，X_imag ，是参考信号的实部和虚部，而Y_real ，Y_imag 是采集信号的实部和虚部。将上文中对于互功率谱的公式展开计算：
$$\begin{gathered} R_{xy}=X(f)Y^{\ast }(f)=(X_{real}+j{X}_{imag})(Y_{real}-j{Y}_{imag})=> \\ {X}_{real}{Y}_{real}-X_{real}j{Y}_{imag}+j{X}_{imag}{Y}_{real}-j{X}_{imag}j{Y}_{imag}=> \\ (X_{real}{Y}_{real}+X_{imag}{Y}_{imag})+(j{X}_{imag}{Y}_{real}-X_{real}j{Y}_{imag}) \end{gathered}$$
即对应代码中

        const double cross_real = X_real * Y_real + X_imag * Y_imag;
        const double cross_imag = X_imag * Y_real - X_real * Y_imag;

4.互功率谱进行归一化防溢出

为了只保留相位信息，去除幅频信息，这个时候计算出模值，再除以模值，即去除幅频信息，即幅频一直为1。

        // 相位归一化（避免除零）
        const double mag = std::hypot(cross_real, cross_imag);
        const double eps = 1e-10;  // 防止除以零

        R[i][0] = cross_real / (mag + eps);
        R[i][1] = cross_imag / (mag + eps);

其中防止模值为0，加上eps

4.对互功率谱进行逆傅里叶变换

通过逆傅里叶变换得到互功率谱的时域图像

// 创建IFFT计划（注意：输出是实数数组）
    fftw_plan planIFFT = fftw_plan_dft_c2r_1d(N, R, cross_corr.data(), FFTW_ESTIMATE);

    // 执行逆FFT
    fftw_execute(planIFFT);

    // 归一化（FFTW逆变换不自动缩放）
    for (int i = 0; i < N; ++i) {
        cross_corr[i] /= N;
    }

找到最大相关峰

此时时域峰值最大的地方，表示相位差最大，即延迟最大。

// 寻找最大相关峰
    int max_idx = 0;
    double max_val = cross_corr[0];
    for (int i = 1; i < N; ++i) {
        if (cross_corr[i] > max_val) {
            max_val = cross_corr[i];
            max_idx = i;
        }
    }

    // 处理负时延：如果峰值在后半部分，返回负索引
    if (max_idx > N/2) {
        max_idx -= N;
    }

二、整体代码与效果

注意：这里的测试文件是两个相同的文件，而可以先将其中一个文件读取掉一部分，以做出两个文件在计算时有明显的延迟，其中最好别超过gcc-phat计算量。

#include "p2p_delay_score.h"
#include <cmath>
#include <complex>
#include <algorithm>
#include <valarray>
#include <numeric>
#include <fftw3.h>

#define GCC_SIZE 800 * 2

int gcc_phat(const std::vector<double>& x, const std::vector<double>& y) {
    const int N = x.size();
    std::vector<double> cross_corr(N, 0.0);

    // 正确分配FFT数组（N点实数FFT生成N/2+1个复数）
    const int N_fft = N / 2 + 1;
    fftw_complex *X = fftw_alloc_complex(N_fft);
    fftw_complex *Y = fftw_alloc_complex(N_fft);
    fftw_complex *R = fftw_alloc_complex(N_fft);

    // 创建FFT计划（注意：输入数组应为非const）
    fftw_plan planX = fftw_plan_dft_r2c_1d(N, const_cast<double*>(x.data()), X, FFTW_ESTIMATE);
    fftw_plan planY = fftw_plan_dft_r2c_1d(N, const_cast<double*>(y.data()), Y, FFTW_ESTIMATE);

    // 执行FFT
    fftw_execute(planX);
    fftw_execute(planY);

    // PHAT加权处理
    for (int i = 0; i < N_fft; ++i) {
        // 计算互功率谱
        const double X_real = X[i][0];
        const double X_imag = X[i][1];
        const double Y_real = Y[i][0];
        const double Y_imag = Y[i][1];

        // X * conj(Y)
        const double cross_real = X_real * Y_real + X_imag * Y_imag;
        const double cross_imag = X_imag * Y_real - X_real * Y_imag;

        // 相位归一化（避免除零）
        const double mag = std::hypot(cross_real, cross_imag);
        const double eps = 1e-10;  // 防止除以零

        R[i][0] = cross_real / (mag + eps);
        R[i][1] = cross_imag / (mag + eps);
    }

    // 创建IFFT计划（注意：输出是实数数组）
    fftw_plan planIFFT = fftw_plan_dft_c2r_1d(N, R, cross_corr.data(), FFTW_ESTIMATE);

    // 执行逆FFT
    fftw_execute(planIFFT);

    // 归一化（FFTW逆变换不自动缩放）
    for (int i = 0; i < N; ++i) {
        cross_corr[i] /= N;
    }

    // 寻找最大相关峰
    int max_idx = 0;
    double max_val = cross_corr[0];
    for (int i = 1; i < N; ++i) {
        if (cross_corr[i] > max_val) {
            max_val = cross_corr[i];
            max_idx = i;
        }
    }

    // 处理负时延：如果峰值在后半部分，返回负索引
    if (max_idx > N/2) {
        max_idx -= N;
    }

    // 清理资源
    fftw_destroy_plan(planX);
    fftw_destroy_plan(planY);
    fftw_destroy_plan(planIFFT);
    fftw_free(X);
    fftw_free(Y);
    fftw_free(R);

    return max_idx;
}

bool p2p_delay()
{
    FILE* pRefFile = fopen("./ref.pcm", "rb");
    FILE* pNearFile = fopen("./near.pcm", "rb");

    FILE* pFile = fopen("./recovery.pcm", "wb+");

    if(pRefFile == nullptr || pNearFile == nullptr)
    {
        printf("fail to open file\n");
        return false;
    }

    std::vector<int16_t> vecRef(GCC_SIZE, 0),vecNear(GCC_SIZE, 0);
    std::vector<double> doubleRef(GCC_SIZE, 0.0), doubleNear(GCC_SIZE, 0.0);

    fread(vecNear.data(), 2, 480, pNearFile);
    int index = 0;
    int total = 0;
    double delay = 0.0;
    while(1)
    {
        delay = 0.0f;
        size_t frames = fread(vecRef.data(), 2, GCC_SIZE, pRefFile);
        if(frames != GCC_SIZE)
            break;

        for (size_t i = 0; i < vecRef.size(); ++i) {
            doubleRef[i] = static_cast<double>(vecRef[i]);
        }

        frames = fread(vecNear.data(), 2, GCC_SIZE, pNearFile);
        if(frames != GCC_SIZE)
            break;

        for (size_t i = 0; i < vecNear.size(); ++i) {
            doubleNear[i] = static_cast<double>(vecNear[i]);
        }
        index = gcc_phat(doubleRef, doubleNear);
        delay= index / double(16000);
        printf("delay:%f \n", delay);
    }
    fclose(pRefFile);
    fclose(pNearFile);

    return true;
}

运算 结果如下：

---

总结

从结果上来看其实会偶尔存在误差，这是由于fft的精度问题，在实际使用中受到环境噪声影响可能误差会更大。其实这种方法最好结合频带分割，对于低频、中频、高频给予不同的权重，而这种权重可以根据能量分布来进行计算。本文中将不过多展示了。

反正收藏也不会看，不如帮忙点个赞吧！