Issues in getting FFT to work on Arduino using my own code

Hi,

I'm working on a project which requires the use of FFT to analyze incoming frequency from the analog-read microphone. I want the magnitude of a specific frequency from a bunch of different frequencies that three 50ohm speakers are emitting, namely 1Khz, 3Khz and 5Khz. I tried using the arduinoFFT library but it gives me major peak and not the frequency of interest. Below is the code I came up with. Could you tell what mistake I'm making here? Thanks.

When I run this code, I run into one major issue, namely, the FFT just doesn't work, forget about getting the magnitude. I'm a DSP engineer. I think the math is right but if you could also cross check it, I would be thankful.


#define SAMPLES 128
#define SAMPLING_FREQUENCY 1000

unsigned int sampling_period_us;
unsigned long microseconds;

struct cmpx {
    float real;
    float imag;
};
typedef struct cmpx COMPLEX;

COMPLEX v[SAMPLES];
COMPLEX w[SAMPLES];

float frequency_of_interest = 1000.0;

void setup() {
  Serial.begin(9600);
  sampling_period_us = round(1000000*(1.0/SAMPLING_FREQUENCY));

  for(int i = 0; i < SAMPLES; i++) {
    float angle = -2 * PI * i / SAMPLES;
    w[i].real = cos(angle);
    w[i].imag = sin(angle);
  }
}

void loop() {
  for(int i=0; i<SAMPLES; i++) {
    microseconds = micros();
   
    v[i].real = analogRead(A0);
    v[i].imag = 0;

    //Serial.print("Analog read: "); // Check analog read value
    //Serial.println(v[i].real);

    while(micros() < (microseconds + sampling_period_us)){}
  }

  fft(v, SAMPLES, w);

  for(int i=0; i<SAMPLES; i++) { // Check FFT output
    Serial.print("FFT output: ");
    Serial.println(sqrt(pow(v[i].real, 2) + pow(v[i].imag, 2)));
  }

  int index = round(frequency_of_interest * SAMPLES / SAMPLING_FREQUENCY);
  float magnitude = sqrt(pow(v[index].real, 2) + pow(v[index].imag, 2));

  Serial.print("Magnitude of frequency of interest: "); // Print magnitude of frequency of interest
  Serial.println(magnitude);

  delay(1000);
}




void fft(COMPLEX *Y, int M, COMPLEX *w)  {
  COMPLEX temp1,temp2;            //temporary storage variables
  int i,j,k;                      //loop counter variables
  int upper_leg, lower_leg;       //index of upper/lower butterfly leg
  int leg_diff;                   //difference between upper/lower leg
  int num_stages=0;               //number of FFT stages, or iterations
  int index, step;                //index and step between twiddle factor
  i=1;                            //log(base 2) of # of points = # of stages
  do
  {
    num_stages+=1;
    i=i*2;
  } while (i!=M);

  leg_diff=M/2;                 //starting difference between upper & lower legs
  step=2;                     //step between values in twiddle.h              
  for (i=0;i<num_stages;i++)      //for M-point FFT                 
  {
    index=0;
    for (j=0;j<leg_diff;j++)
    {
      for (upper_leg=j;upper_leg<M;upper_leg+=(2*leg_diff))
      {
        lower_leg=upper_leg+leg_diff;
        temp1.real=(Y[upper_leg]).real + (Y[lower_leg]).real;
        temp1.imag=(Y[upper_leg]).imag + (Y[lower_leg]).imag;
        temp2.real=(Y[upper_leg]).real - (Y[lower_leg]).real;
        temp2.imag=(Y[upper_leg]).imag - (Y[lower_leg]).imag;
        (Y[lower_leg]).real=temp2.real*(w[index]).real-temp2.imag*(w[index]).imag;
        (Y[lower_leg]).imag=temp2.real*(w[index]).imag+temp2.imag*(w[index]).real;
        (Y[upper_leg]).real=temp1.real;
        (Y[upper_leg]).imag=temp1.imag;
      }
      index+=step;
    }
    leg_diff=leg_diff/2;
    step*=2;
  }
  j=0;
  for (i=1;i<(M-1);i++)           //bit reversal for resequencing data*/
  {
    k=M/2;
    while (k<=j)
    {
      j=j-k;
      k=k/2;
    }
    j=j+k;
    if (i<j)
    {
      temp1.real=(Y[j]).real;
      temp1.imag=(Y[j]).imag;
      (Y[j]).real=(Y[i]).real;
      (Y[j]).imag=(Y[i]).imag;
      (Y[i]).real=temp1.real;
      (Y[i]).imag=temp1.imag;
    }
  }
  return;
}

No need to guess whether the math is right. Feed the FFT with a data array containing a calculated sine wave and check whether the power spectrum is what you expect.

Note: this is not the correct way to time your samples, as it does not correctly handle counter overflow (after about 70 minutes of run time):

   while(micros() < (microseconds + sampling_period_us)){}

This does:

   while(micros() - microseconds < sampling_period_us){}

You can also cross check using the FFT from Matlab or similar tool.

What is your arduino board?

Arduino Mega 2560

I did that using a sine wave lookup table, it matches the theoretical data. It just doesn't run on the Arduino correctly.

I'll make the changes. Thanks.

According to Nyquist theorem, the sampling frequency must be at least a two times more than frequency_of_interest.
With your 1000 Hz sampling the maximum freq you could measure is 500 Hz.
For 3Khz and 5Khz you need a sampling with 6 and 10kHz respectively

Please explain what that means, and give an example.

Never mind. I ran this on an Uno, after reducing the array dimensions to fit, and the result is clearly nonsense. I suspect that the math coding is not correct.

The ArduinoFFT library works as expected.

#define SAMPLES 64
#define SAMPLING_FREQUENCY 1000

unsigned int sampling_period_us;
unsigned long microseconds;

struct cmpx {
    float real;
    float imag;
};
typedef struct cmpx COMPLEX;

COMPLEX v[SAMPLES];
COMPLEX w[SAMPLES];

float frequency_of_interest = 1000.0;

void setup() {
  Serial.begin(115200);
  sampling_period_us = round(1000000*(1.0/SAMPLING_FREQUENCY));

  for(int i = 0; i < SAMPLES; i++) {
    float angle = -2 * PI * i / SAMPLES;
    w[i].real = cos(angle);
    w[i].imag = sin(angle);
  }
}

void loop() {
  for(int i=0; i<SAMPLES; i++) {
//    microseconds = micros();
   
    v[i].real = 100.0*sin(2*PI*i*15./SAMPLES); //15 Hz
    v[i].imag = 0;

    Serial.print("Analog read: "); // Check analog read value
    Serial.println(v[i].real);

//    while(micros() < (microseconds + sampling_period_us)){}
  }

  fft(v, SAMPLES, w);

  for(int i=0; i<SAMPLES; i++) { // Check FFT output
    Serial.print("FFT output: ");
    Serial.println(sqrt(pow(v[i].real, 2) + pow(v[i].imag, 2)));
  }

  int index = round(frequency_of_interest * SAMPLES / SAMPLING_FREQUENCY);
  float magnitude = sqrt(pow(v[index].real, 2) + pow(v[index].imag, 2));

  Serial.print("Magnitude of frequency of interest: "); // Print magnitude of frequency of interest
  Serial.println(magnitude);

  while(1);
}




void fft(COMPLEX *Y, int M, COMPLEX *w)  {
  COMPLEX temp1,temp2;            //temporary storage variables
  int i,j,k;                      //loop counter variables
  int upper_leg, lower_leg;       //index of upper/lower butterfly leg
  int leg_diff;                   //difference between upper/lower leg
  int num_stages=0;               //number of FFT stages, or iterations
  int index, step;                //index and step between twiddle factor
  i=1;                            //log(base 2) of # of points = # of stages
  do
  {
    num_stages+=1;
    i=i*2;
  } while (i!=M);

  leg_diff=M/2;                 //starting difference between upper & lower legs
  step=2;                     //step between values in twiddle.h              
  for (i=0;i<num_stages;i++)      //for M-point FFT                 
  {
    index=0;
    for (j=0;j<leg_diff;j++)
    {
      for (upper_leg=j;upper_leg<M;upper_leg+=(2*leg_diff))
      {
        lower_leg=upper_leg+leg_diff;
        temp1.real=(Y[upper_leg]).real + (Y[lower_leg]).real;
        temp1.imag=(Y[upper_leg]).imag + (Y[lower_leg]).imag;
        temp2.real=(Y[upper_leg]).real - (Y[lower_leg]).real;
        temp2.imag=(Y[upper_leg]).imag - (Y[lower_leg]).imag;
        (Y[lower_leg]).real=temp2.real*(w[index]).real-temp2.imag*(w[index]).imag;
        (Y[lower_leg]).imag=temp2.real*(w[index]).imag+temp2.imag*(w[index]).real;
        (Y[upper_leg]).real=temp1.real;
        (Y[upper_leg]).imag=temp1.imag;
      }
      index+=step;
    }
    leg_diff=leg_diff/2;
    step*=2;
  }
  j=0;
  for (i=1;i<(M-1);i++)           //bit reversal for resequencing data*/
  {
    k=M/2;
    while (k<=j)
    {
      j=j-k;
      k=k/2;
    }
    j=j+k;
    if (i<j)
    {
      temp1.real=(Y[j]).real;
      temp1.imag=(Y[j]).imag;
      (Y[j]).real=(Y[i]).real;
      (Y[j]).imag=(Y[i]).imag;
      (Y[i]).real=temp1.real;
      (Y[i]).imag=temp1.imag;
    }
  }
  return;
}

Thanks, I'll look into the math.

This topic was automatically closed 180 days after the last reply. New replies are no longer allowed.