[linux-audio-dev] Using FFT to figure out frequencies of partials?

[Mario Lang]

Hi.

I am trying to finally get a grip on what a FFT really does by the
learning by doing approach.  I ripped some code from qjacktuner (thanks to
Karsten Wiese) and rewrote it such that it uses fftw3 (which I'd
prefer as a backend for now).  However, the code executes, but returns
very strange results, and so I am asking the DSP coders of you
if you perhaps can see the probably quite obvious error I am making?

When I pluck the A string on my guitar, I receive output like this:

Sample rate: 48000
0 = {freq=143.807441, amp=8.776460}
1 = {freq=162.811147, amp=9.086173}
2 = {freq=181.814852, amp=9.604586}
3 = {freq=200.818557, amp=9.631306}
0 = {freq=148.405443, amp=9.330607}
1 = {freq=271.871957, amp=nan}
2 = {freq=148.405443, amp=49.469994}
3 = {freq=148.405443, amp=45.180032}
0 = {freq=19.003705, amp=47.543385}
2 = {freq=19.003705, amp=66.636058}
3 = {freq=394.003705, amp=63.710202}
0 = {freq=19.003705, amp=19.350767}
2 = {freq=19.003705, amp=28.120214}
3 = {freq=394.003705, amp=56.406310}
0 = {freq=19.003705, amp=42.607500}
2 = {freq=19.003705, amp=78.669607}
3 = {freq=241.974062, amp=38.201785}
0 = {freq=19.003705, amp=2.835922}
2 = {freq=241.974062, amp=93.692646}
3 = {freq=394.003705, amp=48.011389}
0 = {freq=19.003705, amp=45.488988}
2 = {freq=19.003705, amp=46.009611}
3 = {freq=394.003705, amp=60.820928}
2 = {freq=19.003705, amp=54.438861}
3 = {freq=394.003705, amp=47.402751}

In fact, it seems no matter what sound input I give it, it ends up reporting
19.003705 as the loudest frequency.

This is obviously wrong.  Insight appreciated :-)

/* foolowing code is only partial, so it doesn't compile, I can provide
 * you with the full source of the program if you want.
 * Variable "rate" is the current sampling rate (48000 in my case).
 */

#include <fftw3.h>

#define M_PI 3.14159265358979323846
#define MAX_FFT_LENGTH 16384
static float sampleFIFO[MAX_FFT_LENGTH];
static float *sample = NULL;
static float lastPhase[MAX_FFT_LENGTH/2+1];
static int fftSize;
static double *fftIn;
static fftw_complex *fftOut;
static fftw_plan fftPlan;

void
fftInit (int size)
{
  fftIn = fftw_malloc(sizeof(double) * size);
  fftOut = fftw_malloc(sizeof(fftw_complex) * (size/2+1));
  fftSize = size;
  fftPlan = fftw_plan_dft_r2c_1d(fftSize, fftIn, fftOut, FFTW_MEASURE);

  memset(sampleFIFO, 0, MAX_FFT_LENGTH*sizeof(float));
  memset(lastPhase, 0, (MAX_FFT_LENGTH*sizeof(float))/2);
}

void
fftMeasure (int nframes, int osamp, float *indata)
{
  double freqPerBin, expct;
  long i,k, qpd, inFifoLatency, stepSize;

  stepSize = fftSize/osamp;
  freqPerBin = rate/(double)fftSize;
  expct = 2.*M_PI*(double)stepSize/(double)fftSize;
  inFifoLatency = fftSize-stepSize;

  if (!sample) sample = sampleFIFO + inFifoLatency;

  for (i=0; i<nframes; i++) {
    *sample++ = indata[i];

    if (sample-sampleFIFO >= fftSize) {
      sample = sampleFIFO + inFifoLatency;

      for (k=0; k<fftSize; k++) {
        double window = -.5*cos(2.*M_PI*(double)k/(double)fftSize)+.5;
        fftIn[k] = sampleFIFO[k] * window;
      }
      fftw_execute(fftPlan);

      for (k=0; k<=fftSize/2; k++) {
	double magn, phase, tmp, real, imag;

        real = fftOut[k][0];
        imag = fftOut[k][1];

        /* compute magnitude and phase */
        magn = 2.*sqrt(real*real + imag*imag);
        phase = atan2(imag,real);

        /* compute phase difference */
        tmp = phase - lastPhase[k];
        lastPhase[k] = phase;

        /* subtract expected phase difference */
        tmp -= (double)k*expct;

        /* map delta phase into +/- Pi interval */
        qpd = (long)(tmp / M_PI);
        if (qpd >= 0) qpd += qpd&1;
        else qpd -= qpd&1;
        tmp -= M_PI*(double)qpd;

        /* get deviation from bin frequency from the +/- Pi interval */
        tmp = osamp*tmp/(2.*M_PI);

        /* compute the k-th partials' true frequency */
        tmp = (double)k*freqPerBin + tmp*freqPerBin;

	if (tmp > 0.0) {
	  printf("%d = {freq=%f, amp=%f}\n", k, tmp, magn);
	}
      }

      /* move input FIFO */
      for (k=0; k<inFifoLatency; k++) sampleFIFO[k] = sampleFIFO[k+stepSize];
    }
  }
}

void
fftFree ()
{
  fftw_destroy_plan(fftPlan);
  fftw_free(fftIn); fftw_free(fftOut);
}

-- 
Thanks,
  Mario