Go Down

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

#### Judd_Foster ##### Jan 26, 2014, 01:46 pm
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 8float 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 = input + input;  X = input - input;  X = input + input;  X = input - input;  X = input + input;  X = input - input;  X = input + input;  X = input - input;  // butterfly 4  Y = X + X;  Y = X + X;  Y = X - X;  Y = X - X;  Y = X + X;  Y = X + X;  Y = X - X;  Y = X - X;  // final butterfly 8  Z = Y + Y;  Z = Y + (Y * (sina + cosa));  Z = Y + Y;  Z = Y;  Z = Y - Y;  Z = Y - (Y * (sina + cosa));  Z = Y - Y;  Z = Y;    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 #1
##### Jan 26, 2014, 04:33 pm
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 #2
##### Jan 26, 2014, 04:49 pm
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 = 10.00        // input0.711.000.71-0.00-0.71-1.00-0.710.00         // output6.260.000.000.000.000.000.00K = 20.001.00-0.00-1.000.001.00-0.00-1.000.000.000.000.004.000.000.000.00K = 30.000.71-1.000.71-0.00-0.711.00-0.710.000.580.000.000.002.250.000.00`

Also, cosine waveforms produce wacky results:
Code: [Select]
`FFT Values Return:K = 11.000.71-0.00-0.71-1.00-0.710.000.710.001.410.000.590.005.420.000.59K = 21.00-0.00-1.000.001.00-0.00-1.000.000.000.000.000.004.000.004.000.00K = 31.00-0.710.000.71-1.000.71-0.00-0.710.002.590.003.410.000.000.003.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 pmLast 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 = Y + (Y * (sina + cosa));

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

#### Judd_Foster #4
##### Jan 26, 2014, 08:33 pm
What this line
Code: [Select]
`Z = Y + (Y * (sina + cosa));` is supposed to do is use sina and cosa to form
Code: [Select]
`e^2*PI*k*n/N`which is the complex exponential required for FFT. David Dorran said that the e^... value is the sum of sin (Wt)+ cos(Wt):
I got my equations from this video:

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 #5
##### Jan 26, 2014, 09:58 pm
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 #6
##### Jan 27, 2014, 12:03 am
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 = (cosine * inputReal[odd]) - (sine * inputImag[odd])  ti = (wr * fi[j]) + (wi * fr[j]);  // temp = (cosine * inputImag[odd]) + (sine * inputReal[odd])  qr = fr[i];                        // temp = inputReal[even]  qi = fi[i];                        // temp = inputImag[even]  if(shift) {    qr >>= 1;    qi >>= 1;  }  fr[j] = qr - tr;                   // outputReal[odd] = temp - temp  fi[j] = qi - ti;                   // outputImag[odd] = temp - temp  fr[i] = qr + tr;                   // outputReal[even] = temp + temp  fi[i] = qi + ti;                   // outputImag[even] = temp + temp}`

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

Thanks
Judd

#### jremington #7
##### Jan 27, 2014, 01:21 am
That looks better.
tr, ti, qr and qi are the real and imaginary parts of two temporary variables in the original code.

#### Judd_Foster #8
##### Jan 27, 2014, 01:56 am
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?

Judd

#### jremington #9
##### Jan 27, 2014, 02:01 am
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 #10
##### Jan 27, 2014, 02:37 am
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 #11
##### Jan 27, 2014, 02:57 am
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 #12
##### Jan 28, 2014, 12:14 am
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.28float 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 = (input * cosa);  C = (input * cosa);  C = (input * cosa);  C = (input * cosa);  C = (input * cosa);  C = (input * cosa);  C = (input * cosa);  C = (input * cosa);  S = (input * sina);  S = (input * sina);  S = (input * sina);  S = (input * sina);  S = (input * sina);  S = (input * sina);  S = (input * sina);  S = (input * sina);  X = C + S;  X = C - S;  X = C + S;  X = C - S;  X = C + S;  X = C - S;  X = C + S;  X = C - S;  // butterfly 4  Y = (X * cosa) + (X * sina);  Y = (X * cosa) + (X * sina);  Y = (X * cosa) - (X * sina);  Y = (X * cosa) - (X * sina);  Y = (X * cosa) + (X * sina);  Y = (X * cosa) + (X * sina);  Y = (X * cosa) - (X * sina);  Y = (X * cosa) - (X * sina);  // final butterfly 8  output = (Y * cosa) + (Y * sina);  output = (Y * cosa) + (Y * sina);  output = (Y * cosa) + (Y * sina);  output = (Y * cosa) + (Y * sina);  output = (Y * cosa) - (Y * sina);  output = (Y * cosa) - (Y * sina);  output = (Y * cosa) - (Y * sina);  output = (Y * cosa) - (Y * sina);}`
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.001.000.000.000.00-1.00-1.00FFT Values Return:0.00-0.000.00-1.420.00-1.42-0.00-0.00FFT Values Return:0.00-1.00-1.00-0.000.00-0.001.00-1.00FFT Values Return:0.00-0.00-0.000.000.000.000.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 #13
##### Jan 28, 2014, 01:45 am
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 #14
##### Jan 28, 2014, 03:23 am
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