Go Down

Topic: Arduino FFT code not expected results... (Read 7816 times) previous topic - next topic

Judd_Foster

Hello everyone,
I'm trying to make a sketch for Arduino to do the FFT of an 8 sample signal. I looked up on the internet on how to do FFT and came across dozens of webpages explaining what a "butterfly" is. I tried to code up a sketch to do FFT, but I'm not quite sure if it is working or not. When I put certain sinewaves directly in the input[] value, FFT returns expected results, sometimes...
It took me a while just to grasp the basic concept of FFT over DFT, so I don't even know if my math here is right! But here is my code:
Code: [Select]
#define FFTSIZE 8

float input[FFTSIZE];
float output[FFTSIZE];

const float sina[FFTSIZE] = {0, 0.71, 1, 0.71, 0, -0.71, -1, -0.71};
const float cosa[FFTSIZE] = {1, 0.71, 0, -0.71, -1, -0.71, 0, 0.71};

void setup() {
  Serial.begin(9600);
  for(int counter = 0; counter < FFTSIZE; counter++) {
    input[counter] = sin(2 * PI * 1 * counter / FFTSIZE);
    Serial.println(input[counter]);
  }
  Serial.println("FFT Values Return:");
}

void loop() {
  FFT();
  Serial.println();
  for(int counter = 0; counter < FFTSIZE; counter++) {
    Serial.println(output[counter]);
  }
  while(1);
}

void FFT() {
  float X[FFTSIZE];
  float Y[FFTSIZE];
  float Z[FFTSIZE];
  // buterfly 2
  X[0] = input[0] + input[4];
  X[1] = input[0] - input[4];
  X[2] = input[1] + input[5];
  X[3] = input[1] - input[5];
  X[4] = input[2] + input[6];
  X[5] = input[2] - input[6];
  X[6] = input[3] + input[7];
  X[7] = input[3] - input[7];
  // butterfly 4
  Y[0] = X[0] + X[2];
  Y[1] = X[1] + X[3];
  Y[2] = X[0] - X[2];
  Y[3] = X[1] - X[3];
  Y[4] = X[4] + X[6];
  Y[5] = X[5] + X[7];
  Y[6] = X[4] - X[6];
  Y[7] = X[5] - X[7];
  // final butterfly 8
  Z[0] = Y[0] + Y[4];
  Z[1] = Y[1] + (Y[5] * (sina[1] + cosa[1]));
  Z[2] = Y[2] + Y[6];
  Z[3] = Y[3];
  Z[4] = Y[0] - Y[4];
  Z[5] = Y[1] - (Y[5] * (sina[1] + cosa[1]));
  Z[6] = Y[2] - Y[6];
  Z[7] = Y[3];
 
  for(int counter = 0; counter < FFTSIZE; counter++) {
    if(Z[counter] < 0) Z[counter] = 0;
    output[counter] = Z[counter];
  }
}

Any help would be great, thanks.
Judd

Grumpy_Mike

Quote
When I put certain sinewaves directly in the input[] value, FFT returns expected results,

So what is the problem?

Quote
sometimes..

That will be due to the input waveform not being synchronized with a whole number of the sampling rate.

Judd_Foster

There is a number of things wrong with it. Here, I will show you the results of sinewaves in the serial terminal:
Code: [Select]

FFT Values Return:
K = 1
0.00        // input
0.71
1.00
0.71
-0.00
-0.71
-1.00
-0.71

0.00         // output
6.26
0.00
0.00
0.00
0.00
0.00
0.00


K = 2
0.00
1.00
-0.00
-1.00
0.00
1.00
-0.00
-1.00

0.00
0.00
0.00
0.00
4.00
0.00
0.00
0.00


K = 3
0.00
0.71
-1.00
0.71
-0.00
-0.71
1.00
-0.71

0.00
0.58
0.00
0.00
0.00
2.25
0.00
0.00


Also, cosine waveforms produce wacky results:
Code: [Select]

FFT Values Return:
K = 1
1.00
0.71
-0.00
-0.71
-1.00
-0.71
0.00
0.71

0.00
1.41
0.00
0.59
0.00
5.42
0.00
0.59


K = 2
1.00
-0.00
-1.00
0.00
1.00
-0.00
-1.00
0.00

0.00
0.00
0.00
0.00
4.00
0.00
4.00
0.00


K = 3
1.00
-0.71
0.00
0.71
-1.00
0.71
-0.00
-0.71

0.00
2.59
0.00
3.41
0.00
0.00
0.00
3.41

Quote
That will be due to the input waveform not being synchronized with a whole number of the sampling rate.

Could you explain the synchronized part a little more? I don't think I have a sampling rate in my sketch, I haven't hooked up any analog sensors yet.
Thanks,
Judd

jremington

#3
Jan 26, 2014, 06:37 pm Last Edit: Jan 26, 2014, 07:28 pm by jremington Reason: 1
I don't think your math is correct. For example, this line does not look right:
Quote
Z[1] = Y[1] + (Y[5] * (sina[1] + cosa[1]));

Compare your approach with that here: http://wiki.openmusiclabs.com/wiki/ArduinoFFT

Judd_Foster

What this line
Code: [Select]
Z[1] = Y[1] + (Y[5] * (sina[1] + cosa[1])); is supposed to do is use sina[1] and cosa[1] to form
Code: [Select]
e^2*PI*k*n/Nwhich is the complex exponential required for FFT. David Dorran said that the e^... value is the sum of sin (Wt)+ cos(Wt):
http://www.youtube.com/watch?v=K_C7htSXORY
I got my equations from this video:
http://www.youtube.com/watch?v=D5ueRUyCP58
Skip to 7:57 to see 8 sample butterfly.

The Arduino FFT library is written in assembler... I have no background in that, but I will keep trying to make sense of the commands.
Thanks,
Judd

jremington

You can't add the sine and cosine terms like that. They form the two parts of a complex number and have to be treated separately as two real numbers (or you must use the complex variable type in the C/C++ programming language).

I looked around and could not find a good tutorial on the FFT for Arduino. There are many good general tutorials, but a reasonable base for comparison is the  Arduino FIX_FFT (fixed point) code. It is simple and fast and written in plain C. You can find the code and a discussion of some examples of its use at this site: http://members.shaw.ca/el.supremo/Arduino/fft/ffttest.htm

This is the heart of the FFT algorithm (the wr and wi are the real cosine and imaginary sine terms, FIX_MPY is a multiply operation). Just look at the symmetry of the operation.
Code: [Select]
for (i=m; i<n; i+=istep) {
    j = i + l;
    tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
    ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
    qr = fr[i];
    qi = fi[i];
    if (shift) {
  qr >>= 1;
  qi >>= 1;
    }
    fr[j] = qr - tr;
    fi[j] = qi - ti;
    fr[i] = qr + tr;
    fi[i] = qi + ti;
}

Judd_Foster

Is what I put in comments basically correct?

Code: [Select]

for(i = m; i < n; i += istep) {
  j = i + l;
  tr = (wr * fr[j]) - (wi * fi[j]);  // temp[0] = (cosine * inputReal[odd]) - (sine * inputImag[odd])
  ti = (wr * fi[j]) + (wi * fr[j]);  // temp[1] = (cosine * inputImag[odd]) + (sine * inputReal[odd])
  qr = fr[i];                        // temp[2] = inputReal[even]
  qi = fi[i];                        // temp[3] = inputImag[even]
  if(shift) {
    qr >>= 1;
    qi >>= 1;
  }
  fr[j] = qr - tr;                   // outputReal[odd] = temp[2] - temp[0]
  fi[j] = qi - ti;                   // outputImag[odd] = temp[3] - temp[1]
  fr[i] = qr + tr;                   // outputReal[even] = temp[2] + temp[0]
  fi[i] = qi + ti;                   // outputImag[even] = temp[3] + temp[1]
}


If so, can I translate this:
Code: [Select]
Z[1] = Y[1] + (Y[5] * (sina[1] + cosa[1]));
into this:
Code: [Select]
Z[1] = (Y[1] * cosa[1]) + (Y[5] * sina[1]);

Thanks
Judd

jremington

That looks better.
tr, ti, qr and qi are the real and imaginary parts of two temporary variables in the original code.

Judd_Foster

OK great! Now I should do this for all other lines of the fft code, since they all use sina[] and cosa[] in some way, I just replaced the values at the beginning to make the code simpler. I would replace those lines with the new lines right?

Thanks for all your help!
Judd

jremington

Yes, it appears that you optimized away some multiply operations by sines and cosines that happened to be +1 or -1. The code is not easy to follow when you do that.

Judd_Foster

Sorry about that...
I have one more question however. Is there any way to calculate the phase shift of the signal? Like in DFT, when you have the seperate sides of the equation dealing with sin and cos. Is that possible here in FFT?
Thanks,
Judd

jremington

Yes. The amplitude of a particular frequency component is sqrt(re*re + im*im), where re and im are the real and imaginary components of the complex data array element at the frequency index of interest The phase angle in radians is atan2(im,re). When calculating amplitudes, don't forget that the FFT operation scales the transform (by a factor of N if I recall correctly).

Judd_Foster

It still doesn't seem to be working. I replaced all the lines of code with the new good ones. Here is the new and full code:
Code: [Select]

#define FFTSIZE 8
#define _PI 6.28

float input[FFTSIZE];
float output[FFTSIZE];

const float sina[FFTSIZE] = {
  0, 0.71, 1, 0.71, 0, -0.71, -1, -0.71
};
const float cosa[FFTSIZE] = {
  1, 0.71, 0, -0.71, -1, -0.71, 0, 0.71
};

void setup() {
  Serial.begin(9600);
}

void loop() {
  for(int count = 1; count <= 4; count++) {
    for(int counter = 0; counter < FFTSIZE; counter++) {
      input[counter] = sin(_PI * count * counter / FFTSIZE);
    }
    Serial.println();
    FFT();
    Serial.println("FFT Values Return:");
    for(int counter = 0; counter < FFTSIZE; counter++) {
      Serial.println(output[counter]);
    }
  }
  while(1);
}

void FFT() {
  float C[FFTSIZE];
  float S[FFTSIZE];
  float X[FFTSIZE];
  float Y[FFTSIZE];
  float Z[FFTSIZE];
  // buterfly 2
  C[0] = (input[0] * cosa[0]);
  C[1] = (input[0] * cosa[4]);
  C[2] = (input[1] * cosa[0]);
  C[3] = (input[1] * cosa[4]);
  C[4] = (input[2] * cosa[0]);
  C[5] = (input[2] * cosa[4]);
  C[6] = (input[3] * cosa[0]);
  C[7] = (input[3] * cosa[4]);
  S[0] = (input[4] * sina[0]);
  S[1] = (input[4] * sina[4]);
  S[2] = (input[5] * sina[0]);
  S[3] = (input[5] * sina[4]);
  S[4] = (input[6] * sina[0]);
  S[5] = (input[6] * sina[4]);
  S[6] = (input[7] * sina[0]);
  S[7] = (input[7] * sina[4]);
  X[0] = C[0] + S[0];
  X[1] = C[1] - S[1];
  X[2] = C[2] + S[2];
  X[3] = C[3] - S[3];
  X[4] = C[4] + S[4];
  X[5] = C[5] - S[5];
  X[6] = C[6] + S[6];
  X[7] = C[7] - S[7];
  // butterfly 4
  Y[0] = (X[0] * cosa[0]) + (X[2] * sina[0]);
  Y[1] = (X[1] * cosa[2]) + (X[3] * sina[2]);
  Y[2] = (X[0] * cosa[0]) - (X[2] * sina[0]);
  Y[3] = (X[1] * cosa[2]) - (X[3] * sina[2]);
  Y[4] = (X[4] * cosa[0]) + (X[6] * sina[0]);
  Y[5] = (X[5] * cosa[2]) + (X[7] * sina[2]);
  Y[6] = (X[4] * cosa[0]) - (X[6] * sina[0]);
  Y[7] = (X[5] * cosa[2]) - (X[7] * sina[2]);
  // final butterfly 8
  output[0] = (Y[0] * cosa[0]) + (Y[4] * sina[0]);
  output[1] = (Y[1] * cosa[1]) + (Y[5] * sina[1]);
  output[2] = (Y[2] * cosa[2]) + (Y[6] * sina[2]);
  output[3] = (Y[3] * cosa[3]) + (Y[7] * sina[3]);
  output[4] = (Y[0] * cosa[0]) - (Y[4] * sina[0]);
  output[5] = (Y[1] * cosa[1]) - (Y[5] * sina[1]);
  output[6] = (Y[2] * cosa[2]) - (Y[6] * sina[2]);
  output[7] = (Y[3] * cosa[3]) - (Y[7] * sina[3]);
}

I get these results from the serial terminal, each time the loop repeats the input sinewave increases by 1 cycle over the 8 samples:
Code: [Select]

FFT Values Return:
0.00
-1.00
1.00
0.00
0.00
0.00
-1.00
-1.00

FFT Values Return:
0.00
-0.00
0.00
-1.42
0.00
-1.42
-0.00
-0.00

FFT Values Return:
0.00
-1.00
-1.00
-0.00
0.00
-0.00
1.00
-1.00

FFT Values Return:
0.00
-0.00
-0.00
0.00
0.00
0.00
0.00
-0.00

What's wrong with this?! Again, it seems as though it is working with some values, and not for others, specifically try #2.
Thanks,
Judd

el_supremo

Why don't you just get your code working with an FFT library first? Then try writing your own.

Pete
Don't send me technical questions via Private Message.

jremington

Quote
What's wrong with this?!

It is very unlikely that anyone on the forum will have the time and/or patience to go through your code line by line and figure out what the problem(s) is/are. You really need to bite the bullet and with the short example that you've got, for one simple input  (I suggest a square wave of just 1's and -1s), go through the code line by line and verify that each line delivers the expected result according to theory. However, keep in mind that your example, even when working properly, will not be very accurate because the sine/cosine table is accurate to only two digits.

Go Up