在C语言中实现Goertzel算法

Question

我正在DSP处理器上实现BFSK频率跳变通信系统。一些论坛成员建议使用Goertzel算法在特定频率上解调频率跳变。我尝试使用C语言实现了Goertzel算法，代码如下：

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,result,real,imag;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q0 = coeff * q1 - q2 + data[i];
        q2 = q1;
        q1 = q0;
    }
    real = (q1 - q2 * cosine);
    imag = (q2 * sine);
    result = sqrtf(real*real + imag*imag);
    return result;
}

我使用一个函数来计算给定数据集在特定频率下的结果，但是我得到的结果不正确。然而，如果我使用相同的数据集并使用MATLAB goertzel()函数计算goertzel结果，则可以完美地获得结果。我使用C语言实现了该算法，并借助于我在互联网上找到的一些在线教程。我只想听听你们的意见是否正确地实现了goertzel算法。

- anshu

---

你需要将你的答案按样本数的一半进行缩放。 - K. Brafford

如果你有一个好的 Matlab 实现和一个糟糕的 C 实现，你可以在 Matlab 版本中添加日志记录来计算循环中间变量的值，然后将它们与 C 版本的等效变量进行比较。 - Oliver Charlesworth

你解决了你的问题吗？ - K. Brafford

嗨，这个解决方案肯定有效，但是我实现中真正的问题是由于我的IDE错误地连接了文件，导致值没有正确返回...感谢你的帮助... - anshu

Answer 1

如果您认为Matlab的实现之所以好，是因为其结果与您数据的DFT或FFT的该频率的结果相匹配，那么很可能是因为Matlab的实现通过缩放因子对结果进行了标准化，就像FFT一样。

更改您的代码以考虑这一点，并查看是否改善了您的结果。请注意，我还更改了函数和结果名称，以反映您的goertzel计算的是幅值，而不是完整的复杂结果，以增加清晰度：

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

    float   scalingFactor = numSamples / 2.0;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q0 = coeff * q1 - q2 + data[i];
        q2 = q1;
        q1 = q0;
    }

    // calculate the real and imaginary results
    // scaling appropriately
    real = (q1 - q2 * cosine) / scalingFactor;
    imag = (q2 * sine) / scalingFactor;

    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}

- K. Brafford

---

嗨，这个解决方案肯定有效，但我的实现真正的问题是由于我的IDE已将文件链接错误，导致值未正确返回......感谢您的帮助.... - anshu

Answer 2

通常情况下，您可以在计算中使用幅值的平方，例如用于音调检测。

一些出色的Goertzel算法示例可在Asterisk PBX DSP代码 Asterisk DSP code (dsp.c) 和spandsp库 SPANDSP DSP Library 中找到。

- Wadester

Answer 3

考虑两个输入示波形：

1）振幅为A，频率为W的正弦波。

2）振幅和频率与1）相同的余弦波。

Goertzel算法应该对这两个提到的输入波形产生相同的结果，但是提供的代码返回不同的值。我认为代码应该按照以下方式进行修订：

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

    float   scalingFactor = numSamples / 2.0;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q2 = q1;
        q1 = q0;
        q0 = coeff * q1 - q2 + data[i];
    }

    // calculate the real and imaginary results
    // scaling appropriately
    real = (q0 - q1 * cosine) / scalingFactor;
    imag = (-q1 * sine) / scalingFactor;

    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}

- Mohammad Hadi

---

我一直看到这个幅度，但我不确定它如何返回频率？是否有一个公式可以使用幅度计算频率（以赫兹为单位）？ - developer01