1

I am using fftw to get spectrum of an Audio Signal. I get Audio Samples in float32 and have also tried other formats using PortAudio. Then I take fft of it using fftw. In all the formats, I have slight drift of frequency peak from actual frequency. My setup is.

  • Signal Generator to supply sinusoidal signal.
  • Sampling rate of 44.1KHz.

I am getting correct/relatively accurate readings till 10KHz. However as I gradually increase the frequency of Generator, I start getting offset. For example above 10KHz, this is what happens.

Actual frequency         Peak on Spectrum
10KHz                    10.5KHz
12KHz                    12.9KHz
14KHz                    15.5KHz
16KHz                    18.2KHz

The fft Code is as like this.

//Take Samples and do Windowing

 for( i=0; i<framesPerBuffer; i++ )
     {          
         samples[i] = in[i];
         fft->fftIn[i] = (0.54-(0.46*cos(scale_fact*i))) * samples[i];
     }

//Zero Padding

  for(i=framesPerBuffer; i<fftSize; i++)
  {
    fft->fftIn[i] = 0.0;
  }

//FFTW Code

{  fftSize = fftSize;
  this->fftSize = fftSize;
    cout << "Plan start " <<  endl;

  outArraySize = fftSize/2+1;
  cout << "fft Processor start \n";
  fftIn = ((double*) fftw_malloc(sizeof(double) * fftSize));
  fftOut = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * outArraySize );
  fftOutAbs = (double*) fftw_malloc(sizeof(double) * outArraySize );

   fftwPlan = fftw_plan_dft_r2c_1d(fftSize, fftIn, fftOut, FFTW_ESTIMATE);


  //  fftwPlan = fftw_plan_dft_r2c_1d(fftSize, fftIn, fftOut, FFTW_MEASURE);
}

//Absolute response

  int n=fftSize/2+1;
  for(int i=0; i < n; i++)
  {
    fftOutAbs[i] = sqrt(fftOut[i][0]*fftOut[i][0] + fftOut[i][1]*fftOut[i][1]);

  }

  for(unsigned int i=0; i < n; i++)
  {   
     mainCureYData[i] = 20.0 * log10(one_over_n * fft->fftOutAbs[i]);
  }

I need some hints as to where/why this problem could be ?

The Hardware setup look fine as the peak shows correct on sndpeek application.

Thanks,

jav321
  • 151
  • 2
  • 11

1 Answers1

0

You don't actually show the code where you calculate the frequency of the peak, but my guess is that the problem is here:

int n=fftSize/2+1;

This should most likely be:

int n=fftSize/2;

See this answer and/or this answer for an explanation of how to determine frequency from FFT output. (TL;DR: f = Fs * i / n)

Community
  • 1
  • 1
Paul R
  • 208,748
  • 37
  • 389
  • 560
  • For that, this code below generates the bins and then this array is plotted with FFT data to get a graph. The highest magnitude bin should be the frequency of signal. Which is what "drifts" away in my case. `code for(int i=0; i < DEFAULT_FFT_SIZE/2+1; i++) { mainCureXData[i] = ((double) i) / (DEFAULT_FFT_SIZE/2+1) * ap->getSampleRate()/2; ` ` – jav321 Jun 12 '15 at 15:21
  • Again it looks like an over-fondness for adding 1 to things. ;-) `DEFAULT_FFT_SIZE/2+1` should most likely be `DEFAULT_FFT_SIZE/2`. See the linked answers in my answer above to get a better explanation of how bin index maps to frequency. Note also that you probably don't want to look at the Nyquist bin. – Paul R Jun 12 '15 at 15:22
  • the number of FFT points is `fftSize`, so I'd expect `f = Fs * i / fftSize`. The OP's variable `n=fftSize/2+1` is correctly used to process the lower half spectrum. It just shouldn't be used for the frequency computation part. – SleuthEye Jun 12 '15 at 15:24
  • @SlethEye: indeed - I was just guessing as the relevant code was not shown. It looks like a different "off by one" error - see comment above. – Paul R Jun 12 '15 at 15:26
  • Ok . Have tried `mainCureXData[i] = (((double) i) * 44100.0) / 8192 ;` where 44.1K is fs and 8192 is fftSize but the issue persists. Intuitively, can +1 cause such behaviour ? – jav321 Jun 12 '15 at 15:49
  • @jav321 with `fftSize=8192`, the "+1" causes a ~0.02% error in the frequency estimate, not the 5-14% error your results show. There probably are other sources of error or distortions. – SleuthEye Jun 12 '15 at 17:51