Go Down

Topic: Is my Arduino really doing FFT? (Read 6616 times) previous topic - next topic

Judd_Foster

Hi,

I just  finished a sketch for my Arduino to do FFT of some music. I think I understand the FFT and the Arduino seems to be doing what its supposed to be doing, but I just want to make sure I have this right. However, my code seems way too simple based on how other people did this and other people's code is way more complicated than this. I'll post a video link and my code here. I just have the output of a mixer going into the Arduino's analog pin via a diode, so it protects my mixer. A speaker is connected to the other channel so that it doesn't interfere with the FFT codework.

Code: [Select]
// FFT Sketch for Arduino

const int FFTSIZE = 128;
const int ROWS = 8;

const int latchPin = 8;
const int clkPin = 11;
const int dataPin = 12;

const int Sinewave[FFTSIZE] = {
  +0,  +1608,  +3212,  +4808,  +6393,  +7962,  +9512, +11039,
  +12539, +14010, +15446, +16846, +18204, +19519, +20787, +22005,
  +23170, +24279, +25329, +26319, +27245, +28105, +28898, +29621,
  +30273, +30852, +31356, +31785, +32137, +32412, +32609, +32728,
  +32767, +32728, +32609, +32412, +32137, +31785, +31356, +30852,
  +30273, +29621, +28898, +28105, +27245, +26319, +25329, +24279,
  +23170, +22005, +20787, +19519, +18204, +16846, +15446, +14010,
  +12539, +11039,  +9512,  +7962,  +6393,  +4808,  +3212,  +1608,
  +0,  -1608,  -3212,  -4808,  -6393,  -7962,  -9512, -11039,
  -12539, -14010, -15446, -16846, -18204, -19519, -20787, -22005,
  -23170, -24279, -25329, -26319, -27245, -28105, -28898, -29621,
  -30273, -30852, -31356, -31785, -32137, -32412, -32609, -32728,
  -32767, -32728, -32609, -32412, -32137, -31785, -31356, -30852,
  -30273, -29621, -28898, -28105, -27245, -26319, -25329, -24279,
  -23170, -22005, -20787, -19519, -18204, -16846, -15446, -14010,
  -12539  -11039,  -9512,  -7962,  -6393,  -4808,  -3212,  -1608
};

int input[FFTSIZE];
int output[FFTSIZE];
long int sum;
int threshold;

int freeRam() {
  extern int __heap_start, *__brkval;
  int v;
  return (int) &v - (__brkval == 0 ? (int) &__heap_start : (int) __brkval);
}

void setup() {
  pinMode(13, OUTPUT);
  pinMode(latchPin, OUTPUT);
  pinMode(clkPin, OUTPUT);
  pinMode(dataPin, OUTPUT);
  calibAudio();
}

void loop() {
  for(byte counter = 0; counter < FFTSIZE; counter++) {
    input[counter] = analogRead(A0) - threshold;
    if(input[counter] < 5 && input[counter] > -5) {
      input[counter] = 0;
    }
  }
  fftExecute(4);
  ledDisplay(sum, 8);
  fftExecute(7);
  ledDisplay(sum, 7);
  fftExecute(10);
  ledDisplay(sum, 6);
  fftExecute(20);
  ledDisplay(sum, 5);
  fftExecute(30);
  ledDisplay(sum, 4);
  fftExecute(40);
  ledDisplay(sum, 3);
  fftExecute(50);
  ledDisplay(sum, 2);
  fftExecute(60);
  ledDisplay(sum, 1);
}

void calibAudio() {
  digitalWrite(13, HIGH);
  delay(2000);
  for(byte counter = 0; counter < 25; counter++) {
    sum += analogRead(A0);
    delay(1);
  }
  sum /= 25;
  threshold = sum;
  delay(2000);
  digitalWrite(13, LOW);
}

int fftExecute(byte SinewaveNumber) {
  sum = 0;
  long int SinewaveCounter = -1;
  for(byte counter = 0; counter < FFTSIZE; counter++) {
    SinewaveCounter = (SinewaveCounter++) * SinewaveNumber;
    if(SinewaveCounter >= FFTSIZE) {
      SinewaveCounter = 0;
    }
    output[counter] = input[counter] * (Sinewave[SinewaveCounter] / 1024);
  }
  for(byte counter = 0; counter < FFTSIZE; counter++) {
    sum += output[counter];
  }
  if(sum < 0) {
    sum = 0;
  }
  sum = map(sum, 0, 63000, 0, 1024);
  if(sum < 9) {
    sum = 0;
  }
}

void ledDisplay(int sum, int row) {
  // row = (map(row, 0, FFTSIZE / 2, 0, ROWS) + 1);
  sum = map(sum, 0, 1024, 0, 8);
  if(row == 0) {
    row = B00000000;
  }
  else if(row == 1) {
    row = B00000001;
  }
  else if(row == 2) {
    row = B00000010;
  }
  else if(row == 3) {
    row = B00000100;
  }
  else if(row == 4) {
    row = B00001000;
  }
  else if(row == 5) {
    row = B00010000;
  }
  else if(row == 6) {
    row = B00100000;
  }
  else if(row == 7) {
    row = B01000000;
  }
  else if(row == 8) {
    row = B10000000;
  }
  if(sum == 0) {
    sum = B00000000;
  }
  else if(sum == 1) {
    sum = B00000001;
  }
  else if(sum == 2) {
    sum = B00000011;
  }
  else if(sum == 3) {
    sum = B00000111;
  }
  else if(sum == 4) {
    sum = B00001111;
  }
  else if(sum == 5) {
    sum = B00011111;
  }
  else if(sum == 6) {
    sum = B00111111;
  }
  else if(sum == 7) {
    sum = B01111111;
  }
  else if(sum >= 8) {
    sum = B11111111;
  }
  digitalWrite(latchPin, LOW);
  shiftOut(dataPin, clkPin, MSBFIRST, row);
  shiftOut(dataPin, clkPin, MSBFIRST, (255 - sum));
  shiftOut(dataPin, clkPin, MSBFIRST, (255 - sum));
  shiftOut(dataPin, clkPin, MSBFIRST, 255);
  digitalWrite(latchPin, HIGH);
}


https://docs.google.com/file/d/0B2zuATFj02CVZE1oQm1uRF9XRE0/edit?usp=drive_web

I thought that my analogRead function was sampling at 10KHz, but the FFT seems to be running really slowly on my Arduino. The FFT should have been able to pick up frequencies up to 5Khz, but it seems that when I turn the LOW section on my EQ entirely off on my mixer, the entire display just shuts off, no matter how high I crank up the MID and HIGH. This leads me to believe that my FFT sketch isn't doing any analysis on frequencies above around 500Hz. What's going on there? I also have another question about why the third light in from the left on my LED matrix never lights up. The LED and circuitry works just fine, its just in the code, when I use fftExecute(eight) <--- (the #8 would have made smiley-cool) all the way up to fftExecute(12) it doesn't light up any LED's on the display.

Thanks,

Judd

el_supremo

That doesn't look like an FFT. Do you have a reference to whatever method is being used in fftExecute?

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

Judd_Foster

http://www.youtube.com/watch?v=B2iUDBZzBpY at 13:00

When he said to multiply the sinewaves, that made sense to me, and when I tried that using hard-coded values for input[] the entire thing worked perfectly. I don't know if its FFT or DFT, but it's doing something!

Judd

Magician

No, it's has nothing to do with FFT.  "for" loop is summing up multiple of passing constant a a number 0 - 127.
Code: [Select]
int fftExecute(byte SinewaveNumber) {
  sum = 0;
  long int SinewaveCounter = -1;
  for(byte counter = 0; counter < FFTSIZE; counter++) {
    SinewaveCounter = (SinewaveCounter++) * SinewaveNumber;
    if(SinewaveCounter >= FFTSIZE) {
      SinewaveCounter = 0;
    }

BTW, presents sinewave lookup table in the code isn't an real indicator of the fft, butterfly is.
Check on a link below. 

el_supremo

Quote
I don't know if its FFT or DFT

It's neither. I looked at the video and your code and you aren't doing what is described in the video, plus you have several kinds of programming errors. For example, looking at this piece of code:
Code: [Select]
  long int SinewaveCounter = -1;
  for(byte counter = 0; counter < FFTSIZE; counter++) {
    SinewaveCounter = (SinewaveCounter++) * SinewaveNumber;
    if(SinewaveCounter >= FFTSIZE) {
      SinewaveCounter = 0;
    }
    output[counter] = input[counter] * (Sinewave[SinewaveCounter] / 1024);
  }

The first time into the loop, SinewaveCounter is initialized to -1. If we assume that SinewaveNumber is 20, then this statement:
Code: [Select]
    SinewaveCounter = (SinewaveCounter++) * SinewaveNumber;
will set SinewaveCounter to -20 (or it might set it to -19, which is a whole 'nuther problem because of the ++). You then use SinewaveCounter to index into the Sinewave array to access element Sinewave[-20]. But a negative array index is illegal so it will pick a random number in memory rather than an element of the sine wave table.
Then you divide that value by 1024 which introduces more problems. The Sinewave table is an array of integers and you are dividing it by an integer so the division is done as an integer division. The result of this division will only be a number in the range -16 to 15.
Then the next time through the loop, SinewaveCounter will be -400 and things get worse. There is no way that this is computing a DFT.
Another problem is that a DFT requires that you compute a sum of products involving sines and separately, a sum of products involving cosines. You can't just do one of them because you need both in order to determine the magnitude of the signal in each frequency bin.

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

Judd_Foster

Yikes! That doesn't sound too good. Tomorrow I'll look into fixing those issues, but what was he describing in the video? If I add a cosine wave value to this (and somehow fix all my issues) will it be more like what he said? Was what he says in the video an algorithm for fft or was it something totally different?

Thanks for all your patience, as you can see, I don't know a lot about Arduino.

Judd

Magician

Quote
I don't know a lot about Arduino.
Nether  C++, you pick-up a topic that may be considered "advanced"..

el_supremo

What he described in the video is the Discrete Fourier Transform. If you are going to implement it, you will need to have a better understanding of the equation which shows how to compute the DFT. At about 9:30 in the video, the equation on the third line (in blue) is the one you need to implement.
As Magician says, you also need to have a better understanding of how to program in C/C++.
The easiest way out would be to try to find a library that does the DFT for you.

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

Judd_Foster

Magician - I had a look at your link before I even started this project, and that's where I got my inspiration from. It seemed a little complicated for me to understand so I tried to make my own sketch. I'll try to find some tutorials and learn more about the limitations and programming techniques of arduino.

El-supremo - I'll have a look at that equation and do some more learning on programming because what I wanted to do was to implement a sound localizor, like magicians but modify it so that would work on the UNO, and then build a voice recognition system around that based on an SD card.

Thanks guys for your help, btw I'm 13 years old so all this math seems a little daunting to me.

Judd

Judd_Foster

Quote
At about 9:30 in the video, the equation on the third line (in blue) is the one you need to implement.

I tried to implement that equation into the Arduino code. I am finding it a little hard to understand "j", which is supposed to be an imaginary number. What I did was remove "j" from the equation, and just do everything else the same. I did some testing, and my results do work for figuring out what sine waves and cosine waves are in the input[] bin. The code isn't fast, but I just need to understand how DFT works to get my project going.

This is the code I used:

Code: [Select]
#define FFTSIZE 128

float input[FFTSIZE];
float output[FFTSIZE];
float sine;
float cosine;
int counter;
int k;

void setup() {
  Serial.begin(9600);
  Serial.println();
  Serial.println(" SIN    COS    SUM");
 
  for(counter = 0; counter < FFTSIZE; counter++) {
    input[counter] = sin((2 * PI * 1 * counter) / FFTSIZE) + sin((2 * PI * 8 * counter) / FFTSIZE);
  } 
 
  for(k = 1; k <= FFTSIZE / 2; k++) {
    cosine = 0;
    sine = 0;
    output[k] = 0;
    for(counter = 0; counter < FFTSIZE; counter++) {
      cosine += input[counter] * cos((2 * PI * k * counter) / FFTSIZE);
    }
    for(counter = 0; counter < FFTSIZE; counter++) {
      sine += input[counter] * sin((2 * PI * k * counter) / FFTSIZE);
    }
    Serial.print(cosine);
    Serial.print(", ");
    Serial.print(sine);
    Serial.print(", ");
    output[k] = sine - cosine;
    Serial.println(output[k]);
  }
}

void loop() {
}


And the results I got out of the Serial monitor:

Code: [Select]
SIN    COS    SUM
-0.00, 64.00, 64.00
0.00, -0.00, -0.00
-0.00, 0.00, 0.00
0.00, -0.00, -0.00
-0.00, -0.00, -0.00
-0.00, 0.00, 0.00
-0.00, -0.00, -0.00
0.00, 64.00, 64.00
-0.00, -0.00, -0.00
0.00, 0.00, -0.00
-0.00, 0.00, 0.00
0.00, -0.00, -0.00
-0.00, 0.00, 0.00
-0.00, 0.00, 0.00
-0.00, -0.00, -0.00
-0.00, -0.00, -0.00
0.00, -0.00, -0.00
-0.00, -0.00, -0.00
0.00, -0.00, -0.00
-0.00, -0.00, 0.00
0.00, 0.00, 0.00
0.00, 0.00, -0.00
0.00, 0.00, -0.00
0.00, -0.00, -0.00
0.00, -0.00, -0.00
-0.00, -0.00, -0.00
0.00, -0.00, -0.00
0.00, -0.00, -0.00
0.00, 0.00, 0.00
-0.00, 0.00, 0.00
-0.00, 0.00, 0.00
-0.00, -0.00, 0.00
0.00, -0.00, -0.00
-0.00, 0.00, 0.00
0.00, 0.00, 0.00
-0.00, 0.00, 0.00
0.00, 0.00, 0.00
0.00, 0.00, -0.00
-0.00, 0.00, 0.00
-0.00, -0.00, 0.00
-0.00, 0.00, 0.00
0.00, -0.00, -0.00
0.00, -0.00, -0.00
-0.00, -0.00, 0.00
0.00, 0.00, -0.00
0.00, 0.00, 0.00
0.00, -0.00, -0.00
-0.00, 0.00, 0.00
-0.00, -0.00, -0.00
0.00, 0.00, 0.00
-0.00, -0.00, -0.00
0.00, 0.00, 0.00
-0.00, -0.00, -0.00
0.00, 0.00, 0.00
0.00, -0.00, -0.00
-0.00, 0.00, 0.00
0.00, -0.00, -0.00
0.00, 0.00, 0.00
0.00, -0.00, -0.00
-0.00, -0.00, -0.00
0.00, -0.00, -0.00
0.00, 0.00, 0.00
0.00, -0.00, -0.00
-0.00, 0.00, 0.00


Please tell me if this is what is supposed to happen or if I have this wrong. I would also like to know how to program "j" into the Arduino, to complete the equation.

Thanks

Judd

el_supremo

That is a lot closer.
A couple of minor things:
- in the loop which does the DFT you start with "for(k=1;" it should be "for(k=0;".
- You print the header " SIN    COS    SUM", but you print the results in the order COS SIN SUM
- the SUM that you print isn't needed. The useful number to print is the magnitude of each bin. Instead of printing sine-cosine you need to print sqrt(sine*sine + cosine*cosine).
- just FYI: the number of samples does not have to be the same as the size of the DFT. The number of samples can be, and usually is, larger than the DFT size.

In the video, he explains at one point that the numbers that you get from the DFT aren't the actual amplitudes of each frequency. You have to divide by (N/2) where N is the number of samples. So, the magnitude of each frequency bin is sqrt(sine*sine + cosine*cosine) but the amplitude of the frequency is sqrt(sine*sine + cosine*cosine)/(N/2). In your sketch the two sines (remember the header is wrong) have a value of 64. Your number of samples is 128 so dividing by 64 gives a result of 1 which is correct because both sine waves have an amplitude of 1.
Fix up your sketch to use "for(k=0;", change the header to " COS    SIN    SUM" and try:
Code: [Select]
    input[counter] = 7*sin((2 * PI * 1 * counter) / FFTSIZE) + 10*cos((2 * PI * 8 * counter) / FFTSIZE);

In the resulting table you will now get one large spike in the SIN column corresponding to the first term (frequency 1) and an even larger spike in the  COS column corresponding to the second term (frequency 8 ). When you use "for(k=0;", you count the lines starting from zero so that frequency 1 will be the second line and frequency 8 will be the ninth line.

Quote
I am finding it a little hard to understand "j"

It is the square root of -1 (and, BTW, in mathematics texts it is represented by "i" instead of "j"). If you're 13 years old, you may have run into the concept of imaginary numbers when solving the roots of a quadratic equation, but it'll be a while before you cover complex numbers which is what the DFT uses. For what you are doing, you don't need to worry about it. The number that you need from the DFT is the amplitude of each frequency and when you compute the magnitude of each frequency bin using sqrt(sine*sine + cosine*cosine), it makes the j disappear.

However, if you want to do sound localization and voice recognition you have a very long way to go.
Quote
The code isn't fast

That is one of the problems you will run into on the UNO. It does not have floating point instructions so calculations involving sin(), cos() and floating point arithmetic are very slow. You aren't going to be able to sample the data fast enough to get a meaningful result without using a lookup table for the sin/cos and other tricks to speed up the process.


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

Judd_Foster

Thank you so much! I have fixed all those silly mistakes you said about before (I woke up early this morning). I didn't know what to do for the sum so what I did was use cosine - sine.
Quote
in the loop which does the DFT you start with "for(k=1;" it should be "for(k=0;".

I purposely didn't use k = 0 because the input would have to be completely silent for it to return a positive value, and if the input was totally silent, none of the other frequency bins would return a positive value. I didn't think it would be useful.
Quote
just FYI: the number of samples does not have to be the same as the size of the DFT. The number of samples can be, and usually is, larger than the DFT size.

The number of samples is actually twice as much as the DFT size, I used this line:
Code: [Select]
for(k = 1; k <= FFTSIZE / 2; k++) {
which should give me only the first 64 frequency bins.

I tried out the code:
Code: [Select]
input[counter] = 7*sin((2 * PI * 1 * counter) / FFTSIZE) + 10*cos((2 * PI * 8 * counter) / FFTSIZE);
Does that just increase the amplitude of the waves?

Quote
Quote

The code isn't fast

That is one of the problems you will run into on the UNO. It does not have floating point instructions so calculations involving sin(), cos() and floating point arithmetic are very slow.

My idea was to hook in an SD card into the Arduino, and so the Arduino could put in the input values to the SD card, and when the input is finished receiving data, the Arduino can pull out the values and do the DFT, so the Arduino wouldn't have to do calculations in real time. Would this be faster while still using an UNO, or should I purchase a better Arduino board?

My next question would be for implementing the analog input for the Arduino. I have used magician's "audio input to Arduino" page at: http://coolarduino.wordpress.com/2012/06/22/audio-input-to-arduino/ and I have moved the codework around in my sketch for it to input audio. I don't know if it works or not, but here is my sketch:
Code: [Select]
#include "TimerOne.h"

#define FFTSIZE 16
#define ROWS 8

float input[FFTSIZE];
float output[FFTSIZE / 2];
float sine;
float cosine;

boolean finished;

static uint8_t  n_sampl;

int16_t temp;

int counter;
int k;
int threshold;
int rowCounter;

const int latchPin = 8;
const int dataPin = 12;
const int clkPin = 11;

int freeRam () {
  extern int __heap_start, *__brkval;
  int v;
  return (int) &v - (__brkval == 0 ? (int) &__heap_start : (int) __brkval);
}

void setup() {
  pinMode(latchPin, OUTPUT);
  pinMode(dataPin, OUTPUT);
  pinMode(clkPin, OUTPUT);
  pinMode(13, OUTPUT);

  Serial.begin(9600);
 
    ADCSRA = 0x87; // freq = 1/128  125 kHz. 13 cycles x 8     usec =  104 usec.
//ADCSRA = 0x86; // freq = 1/64   250 kHz. 13 cycles x 4     usec =   52 usec.
//ADCSRA = 0x85; // freq = 1/32   500 kHz. 13 cycles x 2     usec =   26 usec.
//ADCSRA = 0x84; // freq = 1/16     1 MHz. 13 cycles x 1     usec =   13 usec.
//ADCSRA = 0x83; // freq = 1/8      2 MHz. 13 cycles x 0.5   usec =  6.5 usec.
//ADCSRA = 0x82; // freq = 1/4      4 MHz. 13 cycles x 0.25  usec = 3.25 usec.

  ADMUX = 0x40;
  ADCSRA |= (1 << ADSC);
 
  calibAudio();
 
  Timer1.initialize(150000);
  Timer1.attachInterrupt(audio);

  calibAudio();
}

void loop() {
  audio();
  for(k = 1; k <= FFTSIZE / 2; k++) {
    cosine = 0;
    sine = 0;
    output[k] = 0;
    for(counter = 0; counter < FFTSIZE; counter++) {
      cosine += input[counter] * cos((2 * PI * k * counter) / FFTSIZE);
    }
    for(counter = 0; counter < FFTSIZE; counter++) {
      sine += input[counter] * sin((2 * PI * k * counter) / FFTSIZE);
    }
    output[k] = sine - cosine;
  }

  finished = true;
 
  LEDDisplay();                          // to be filled in later
}

void calibAudio() {
  for(counter = 0; counter < 32; counter++) {
    if(ADCSRA & 0x10) {
      temp = ADCL;
      temp += (ADCH << 8);
      ADCSRA |= (1 << ADSC);
      threshold += temp;
    }
  }
  threshold /= 32;
}

void audio() {
  digitalWrite(13, HIGH);
  if(finished == true) {
    if(ADCSRA & 0x10) {
      temp = ADCL;
      temp += (ADCH << 8);
      ADCSRA |= (1 << ADSC);
      input[n_sampl] = temp - threshold;
    }
    if(++n_sampl >= FFTSIZE) {
      n_sampl = 0;
    }
    finished = false;
  }
  digitalWrite(13, LOW);
}


If there is anything wrong with that please let me know. My end goal is to send all of these variables out into my 8X8 Matrix via 74HC595N shift registers on a breadboard. How would I interface this new codework to my display? Is it the same as my very first code with just new map() commands?

Thank you for all of your help!

Judd

el_supremo

Quote
I purposely didn't use k = 0 because the input would have to be completely silent for it to return a positive value

No. The k=0 term is the DC component of the input signal. With an input signal like this:
Code: [Select]
input[counter] = sin((2 * PI * 1 * counter) / FFTSIZE)
the values stored in the input array will have values which vary between -1 and +1. They vary equally above and below zero so the magnitude of the first bin will be zero. If instead your signal was:
Code: [Select]
input[counter] = 5 + sin((2 * PI * 1 * counter) / FFTSIZE)
the 5 would raise the whole signal so that it is in the range 4 to 6. That would cause a large magnitude in the first frequency bin.

Quote
which should give me only the first 64 frequency bins.

Yeah, I missed that.


Quote
Does that just increase the amplitude of the waves?

Yes. Instead of sin() it has 7*sin() so it will increase the amplitude of the sin wave by a factor of 7. The cosine will be ten times greater and both will be reflected in the output with higher values in the respective bins.
Quote
I tried out the code

so you should have noticed that the amplitudes were higher than before - one in the SIN column and the other in the COS column.

Quote
when the input is finished receiving data, the Arduino can pull out the values and do the DFT

The SD card would only slow things down even more.

Quote
implementing the analog input for the Arduino

I'll let somebody else handle that.

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

Judd_Foster

Quote
No. The k=0 term is the DC component of the input signal. With an input signal like this:
Code: [Select]
input[counter] = sin((2 * PI * 1 * counter) / FFTSIZE)

I don't understand that. Can you put that simpler? The only time I use k is to count up the frequency bins.

Quote
The SD card would only slow things down even more.

How can I make things faster? I'm trying to avoid getting a new Arduino Mega or DUE, but I have tons of UNO's and 1 Micro lying around.

Thanks

Judd

jremington

#14
Nov 11, 2013, 10:12 pm Last Edit: Nov 11, 2013, 10:25 pm by jremington Reason: 1
Use an 8 or 16-bit integer FFT or FHT calculation. Either will be much, much faster. There are several versions available. There are some limitations of the FFT or FHT over the DFT, such as the number of points must be a power of 2.

Edit: here is a collection http://wiki.openmusiclabs.com/wiki/ArduinoFFT

Go Up