如果两个 Goertzel 滤波器的复频率 bin 输出为$a$和$b$，您可以将相位计算为$\alpha = \text{atan2}\left(\text{Im}(a),\ \text{Re}(a)\right)$和$\beta = \text{atan2}\left(\text{Im}(b),\ \text{Re}(b)\right)$， 在哪里$\text{Re}$返回实部和$\text{Im}$返回参数的虚部。大多数编程语言都有两个参数$\text{atan2}$功能。相位差可以计算为$\alpha - \beta = \text{atan2}\left(\text{Im}\left(a\ \text{conj}(b)\right),\ \text{Re}\left(a\ \text{conj}(b)\right)\right)$， 在哪里$\text{conj}$表示取复共轭，为你节省一个$\text{atan2}$以一个复数乘法为代价，对相位进行评估和包装。

还有一个技巧可以计算$a\ \text{conj}(b)$. 乘以 Goertzel 输出$a$通过第二个输入信号和低通滤波结果。在频域中它看起来像这样：

低通滤波隔离对应于$a\ \text{conj}(b)$. 这以低通滤波器为代价为您节省了一个 Goertzel 滤波器。如果你使用这个技巧，你还必须考虑低通滤波器引起的相移。

为了快速$\text{atan2}$请参阅在 FPGA 上计算定点 atan2 的方法。

以下只是我自己的代码中的一个懒惰的复制粘贴，以防您被卡在 goertzel 过滤器的实/图像提取的实现中。否则，去看看奥利的答案。

template<typename Scalar, class Vector>
static std::complex<Scalar> goertzel(const Vector & data, std::size_t size, Scalar omega)
{
    Scalar sine, cosine, coeff, q0(0), q1(0), q2(0);

    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;

    for (Types::fint_t t = 0; t < size; t++)
    {
        q0 = coeff * q1 - q2 + data[t];
        q2 = q1;
        q1 = q0;
    }

    Scalar real = (q1 - q2 * cosine);
    Scalar imag = (q2 * sine);

    return std::complex<Scalar>(real, imag);

}

在哪里...

    angle = std::arg(complex_number); // which is equivalent to atan2(imag, real);
    magnitude = std::abs(complex_number); // again, same as sqrt(imag^2 + real^2);