具有复杂输出或幅度+相位的C/C++ Goertzel算法？答案

【问题标题】：C/C++ Goertzel algorithm with complex output or magnitude+phase?具有复杂输出或幅度+相位的C/C++ Goertzel算法？
【发布时间】：2012-01-12 13:25:59
【问题描述】：

有人知道我在哪里可以获得代码或库来执行具有复杂输出的 Goertzel 算法吗？（或任何其他 1-bin-DFT 算法？）

【问题讨论】：

en.wikipedia.org/wiki/Goertzel_algorithm 呢？
我对信号处理一无所知，但是你看过FFTW吗？
Wikipedia-Implementation 不会创建复杂的输出，只是功率（=magnitude？）
"FFTW 目前没有实现任何通用的剪枝 FFT 算法。"
你在dsp.stackexchange.com上可能会有更好的运气

标签： c++ c algorithm signal-processing

【解决方案1】：

以下是几年前我写的代码。随意使用它，并注明出处。

goertzelfilter.h

/* goertzelfilter.h
*/

#ifndef GOERTZELFILTER_H_
#define GOERTZELFILTER_H_

#include <complex.h>

typedef struct goertzelfilterstruct {
  double coeff ;
  double sine ;
  double cosine ;
} GoertzelFilter;

GoertzelFilter goertzelSetup( double normalizedfreq );

double complex goertzelFilterC( double *sample, int nsamples, GoertzelFilter *g );

#endif

goerzelfilter.c

/* goertzelfilter.c
*/

#include <math.h>
#include <stdlib.h>
#include <complex.h>

#include "goertzelfilter.h"

GoertzelFilter goertzelSetup( double normalizedfreq )
{
  double w = 2*M_PI*normalizedfreq;
  double wr, wi;

  GoertzelFilter g;

  wr = cos(w);
  wi = sin(w);
  g.coeff = 2 * wr;
  g.cosine = wr;
  g.sine = wi;

  return g;
}

double complex goertzelFilterC( double *samples, int nsamples, GoertzelFilter *g )
{
  double sprev = 0.0;
  double sprev2 = 0.0;
  double s, imag, real;
  int n;

  for (n=0; n<nsamples; n++ ) {
    s = samples[n] + g->coeff * sprev - sprev2;
    sprev2 = sprev;
    sprev = s;
  }

  real = sprev*g->cosine - sprev2;
  imag = -sprev*g->sine;

  return real + I*imag;
}

为了测试，你可以试试这个，

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>

#include "goertzelfilter.h"

#define LEN(a) (sizeof(a)/sizeof(a[0]) )

int main() {

  // This will hold our input and output data
  double data[1024] = { 0. };
  double complex filtered[1024] = { 0. };

  // This will hold the filter constants
  GoertzelFilter g = { 0. };

  int n;
  int nwindow = 16;

  double f = 4./LEN(data) ;

  // Generate data with noise
  for ( n = 0 ; n < LEN(data) ; n++ ) {
    data[n] = sin( n * (2.*M_PI) * f ) + 0.5*((float)rand()/RAND_MAX - 0.5);
  }

  // Set up the filter constants, note that we choose a frequency
  g = goertzelSetup( f );

  // Filter the data using a sliding window
  for( n = 0 ; n < LEN(data)-nwindow ; n++ ) {
    filtered[n+nwindow/2]  = goertzelFilterC( &data[n], nwindow, &g )/nwindow;
  }

  // Print the real Valued Data (1st column) and complex valued Goertzel output
  for( n = 0 ; n < LEN(data); n++ ) {
    printf( "%g %g %g\n", data[n], creal(filtered[n]), cimag(filtered[n]) );
  }

}

这是一个显示测试代码的输入和输出的图表：

【讨论】：

使用documenting comments 似乎更有用。
我会考虑一下。函数 goertzelFilterC() 与许多教科书中描述的方式密切相关。常量在 goertzelSetup() 中计算并保存，因此可以以相同的频率重复调用 Goertzel 函数，开销低。当然，另一种用途是只为所有数据调用一次，然后报告相位和功率或幅度。
@grebeard - “记录 cmets”的链接似乎已损坏。
现在有一个“官方”网站，www.doxygen.nl。
@greybeard - 根据您的建议，我在测试代码中添加了 cmets 以解释如何使用 Goertzel 例程，并且添加了显示输出的图表。