14

我正在 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 算法，我只想了解你们的看法。

4

3 回答 3

14

如果您说 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;
}
于 2012-07-20T14:20:34.840 回答
1

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

Goertzel 的一些优秀示例在 Asterisk PBX DSP 代码 Asterisk DSP 代码 (dsp.c) 和 spandsp 库SPANDSP DSP 库中

于 2013-10-08T20:14:14.300 回答
1

考虑两个输入样本波形：

1) 幅度为 A 频率为 W 的正弦波

2）具有相同幅度和频率的余弦波A和W

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;
}
于 2016-12-31T11:11:45.737 回答