我正在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算法，我只想获得你们的观点。
如果你说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;
}
通常你可以在计算中使用幅度的平方，
例如，用于音调检测。
Goertzels的一些优秀示例在Asterisk PBX DSP代码中
星号DSP代码（dsp.c）
并在spandsp库中 SPANDSP DSP库
考虑两个输入样本波形：
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;
}