一阶 Goertzel 算法的另一种方法

Mosi_62

4.54/5 (18投票s)

2013年10月7日

CPOL

4分钟阅读

36439

461

这种方法展示了一种非常容易理解的一阶 Goerzel 算法的实现方式。

下载源代码 - 36.8 KB

一阶 Goertzel 算法

基于傅立叶变换的方程

一阶 Goertzel 算法基于 x[n] 和一个附加信号 h[n] 的卷积，最终，在经过一番复杂的解释后，得到简单的公式

用

和

我首先在一本关于 Matlab 的非常有趣的书籍中发现了这种傅立叶变换的实现。它基本上看起来像这样

a1 = -exp(1i*2*pi*k/N);
y = x(1);
for n = 2:N
y = x(n) - a1*y;
end
y = -a1*y;

好吧，我们现在不分析这个语法。同样的算法在 C 中看起来像这样

for (k = 0; k < N; k++)
{
  c[k].real = y[0].real;
  c[k].imag = y[0].imag;
  w.real = -Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  w.imag = Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  for (n = 1; n <= N; n++)
    c[k] = kdiff(y[n], kprod(c[k], w));
    c[k] = kprod(c[k], w);
    c[k].real = -c[k].real / (double)(N) * 2.0;
    c[k].imag = -c[k].imag / (double)(N) * 2.0;
  }

带有复数运算

public TKomplex kdiff(TKomplex a, TKomplex b)
{
  TKomplex res;
  res.real = a.real - b.real;
  res.imag = a.imag - b.imag;
  return (res);
}
public TKomplex kprod(TKomplex a, TKomplex b)
{
  TKomplex res;
  res.real = a.real * b.real - a.imag * b.imag;
  res.imag = a.real * b.imag + a.imag * b.real;
  return (res);
}

方法的第一步

不知何故，我觉得我好像见过类似的东西

在其他地方…… 有关多项式

y = a1*x + a2*x^2 + a3*x^3 + a4*x^4…

以及根据霍纳法则进行多项式求值的计算。

y = x * (a1 + x*(a2 + x*(a3 + x*(a4 …))))

如果我在一个小循环中实现这一点，它看起来几乎和上面的公式一样。但傅立叶变换看起来不像这个多项式。傅立叶变换的公式是

好的。但有一种方法可以将它变成一个多项式。表达式

是一个长度为 1 的复数向量，绕着单位圆旋转。在复数计算中，我们可以说

并且傅立叶变换可以写成

这可以写成一个多项式。

如果我把它放在一个循环中，它开始于

然后继续像这样

和

直到我们到达 f(0)

理论上

完成这个多项式。

在傅立叶分量的特殊情况下，可以跳过最后一项，因为

在 C 中实现，这看起来像这样

for (k = 0; k < N; k++)
{
  c[k].real = y[N].real;
  c[k].imag = y[N].imag;
  w.real = Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  w.imag = -Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  for (n = N - 1; n > 0; n--)
    c[k] = ksum(y[n], kprod(c[k], w));
    c[k].real = c[k].real / (double)(N) * 2.0;
    c[k].imag = -c[k].imag / (double)(N) * 2.0;
  }

越来越近

这看起来已经很像 Goerzel 算法了。但有三点不同

w 的符号相反。
在内循环中，索引 n 从 N-1 开始，而不是 1。这意味着我们从后往前工作。
在循环中，我们有一个加法而不是减法。

现在最大的问题是：这意味着什么？

如果我在循环中朝另一个方向运行，这对样本 y[n] 的寻址没有任何区别。它们仍然相同。但 w 是一个复数向量。如果我们从 0 开始 n 并将其增加到 N，w 将逆时针旋转。如果我们从 N-1 开始向另一个方向运行，并将 n 递减到 1，w 将顺时针旋转。这会产生一些影响。此外，我们不是从相同的 w 值开始。我们从

而不是

这基本上是 w 的共轭复数值。幸运的是，有了这个起始值，w 将自动顺时针旋转，因为它有一个负的虚部。

最后：由于我们朝另一个方向运行，循环不会在 e⁰ 处结束，我们不能直接跳过最后一项。通过这种修改，算法变为这样

for (k = 0; k < N; k++)
{
  c[k].real = y[0].real;
  c[k].imag = y[0].imag;
  w.real = Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  w.imag = Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  for (n = 1; n <= N; n++)
    c[k] = ksum(y[n], kprod(c[k], w));
    c[k] = kprod(c[k], w);
    c[k].real = c[k].real / (double)(N) * 2.0;
    c[k].imag = -c[k].imag / (double)(N) * 2.0;
  }

更近了

但仍然：Goertzel 算法不使用 e^-jx 的共轭复数。它使用它的负值并减去它而不是加它。因此，最终傅立叶分量的实部变为负值，并且必须在循环结束时取反。现在算法看起来像这样

for (k = 0; k < N; k++)
{
  c[k].real = y[0].real;
  c[k].imag = y[0].imag;
  w.real = -Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  w.imag = Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
  for (n = 1; n <= N; n++)
    c[k] = kdiff(y[n], kprod(c[k], w));
    c[k] = kprod(c[k], w);
    c[k].real = -c[k].real / (double)(N) * 2.0;
    c[k].imag = -c[k].imag / (double)(N) * 2.0;
  }

这就是 Goertzel 算法。从多项式计算推导而来，经过一些小的修改，我们得到了相同的算法。

Goertzel 算法的优点和缺点

Goertzel 算法的主要目标是计算速度。它必须计算比 DFT 的标准实现少得多的正弦和余弦值，因此它的工作速度相当快。但在我关于快速精简 DFT 的文章中，我找到了一种方法来进一步加速 DFT 算法。

我看到了一个可能的缺点：如果我们的样本数量很大，N 也会很大，因此 w 将有一个非常小的角度，必须多次相乘。必须多次处理一个非常小的虚部，与虚部相比，实部很大。这可能会导致这种计算的精度损失。

结论

这种方法展示了一种非常容易理解的一阶 Goerzel 算法的实现方式，并表明该算法或多或少是通过霍纳法则向另一个方向运行的多项式求值。 :-)

但当然：二阶 Goertzel 算法是另一回事。