Difference between revisions of "Fast Fourier Transform Library & Arduino"

From ESE205 Wiki
Jump to: navigation, search
(Process)
Line 52: Line 52:
  
 
Below you can see some classic implementations in other languages, but as you can see the actual implementation is incredibly difficult. You can find an example of a custom FFT at the bottom. Typically one would not code an FFT from scratch unless you had a very specific purpose for it. Many libraries that deal with computational problems regarding any sort of waveform will have an FFT coded in the API that has been vetted by multiple software developers.  
 
Below you can see some classic implementations in other languages, but as you can see the actual implementation is incredibly difficult. You can find an example of a custom FFT at the bottom. Typically one would not code an FFT from scratch unless you had a very specific purpose for it. Many libraries that deal with computational problems regarding any sort of waveform will have an FFT coded in the API that has been vetted by multiple software developers.  
 +
  
 
To implement a general FFT in an Arduino here are the steps:
 
To implement a general FFT in an Arduino here are the steps:
Line 60: Line 61:
  
 
3. run the C++ code given below
 
3. run the C++ code given below
 +
  
 
To really see it in action:
 
To really see it in action:

Revision as of 17:52, 8 December 2018

Overview

Introduction

As the name suggests the Fast Fourier Transform Library enables for the timely computation of a signal's discrete Fourier transform. Instructions on how to download the latest release can be found here. A fast Fourier transform (fFt) would be of interest to any wishing to take a signal or data set from the time domain to the frequency domain.

Materials & Prerequisites

Materials

All items one needs to utilize an FFT with an Arduino are:

    *Computer
    *Arduino
    *USB connector

Prerequisites

Next download latest open-source Arduino Software here.

Process

First, let's begin with a discussion of a general Fourier Transform (FT) and then we will address the Fast Fourier Transform (FFT). The Fourier transform is a mathematical function that decomposes a waveform, which is a function of time, into the frequencies that make it up. The result produced by the Fourier transform is a complex valued function of frequency. The absolute value of the Fourier transform represents the frequency value present in the original function and its complex argument represents the phase offset of the basic sinusoidal in that frequency.

The Fourier transform is also called a generalization of the Fourier series. This term can also be applied to both the frequency domain representation and the mathematical function used. The Fourier transform helps in extending the Fourier series to non-periodic functions, which allows viewing any function as a sum of simple sinusoids. The definition of which is provided below. Note that we use the discrete time definition of the FT, or a Discrete Fourier Transform (DFT), and not the continuous time definition, the main difference between the two being a summation vs an integral. We make this choice because inputs into computers are discrete data points and not continuous so we can not use an integral, and by extension, we can't use the continuous time definition of the FT.

The decomposition of functions into their respective frequencies is a very powerful and useful tool, however, an FT requires massive amounts of computational work. For those of you in computer science, an FT would take O(n^2) time. So although we want to use this tool it is computationally expensive. Or at least this was the case until J.W Cooley and John Tukey came on the scene. They came up with the aptly named Cooley–Tukey FFT algorithm which reduced the time cost on a DFT from O(n^2) to O(nlogn). It recursively breaks down the DFT into smaller DFTs and turns the summation into a dynamic programming problem. So we save on time but we increase our space complexity dramatically, however, due to the continued improvement in transistor technology, modern day computing, for all practical purposes, doesn't care about space complexity until it has to. An example of a field that has to would be computational biology, due to the vast number of different genes in existence that all need to be tested. But that's beside the point, let's continue our discussion of an FFT.

The actual algorithm begins by splitting the matrix into two parts, one with all the even indexed elements, the other part with all the odd indexed elements. We continue spliting the matrices down in the same manner until we can perform a DFT on a manageable matrix. Another way to implement the split, rather than by the even-odd index convention, would be through a reversing the bit value of the array entry's index. You make the choice only on the language you program in; so choose the option that takes up less time to implement for the language.

Coding it will look like:


function cooley_tukey(x)

   N = length(x)
   if (N > 2)
       x_odd = cooley_tukey(x[1:2:N])
       x_even = cooley_tukey(x[2:2:N])
   else
       x_odd = x[1]
       x_even = x[2]
   end
   n = 0:N-1
   half = div(N,2)
   factor = exp.(-2im*pi*n/N)
   return vcat(x_odd .+ x_even .* factor[1:half], x_odd .- x_even .* factor[1:half])

end

Below you can see some classic implementations in other languages, but as you can see the actual implementation is incredibly difficult. You can find an example of a custom FFT at the bottom. Typically one would not code an FFT from scratch unless you had a very specific purpose for it. Many libraries that deal with computational problems regarding any sort of waveform will have an FFT coded in the API that has been vetted by multiple software developers.


To implement a general FFT in an Arduino here are the steps:

1. Get your computer, arduino, usb-B cable ready

2. download the Arduino coding terminal on your computer

3. run the C++ code given below


To really see it in action:

1. Find a digital waveform (an audio signal, voltage signal, anything digital that can be modeled with sinusoids really)

2. Pass that waveform in as the parameter of the FFT

3. Print the output


Implementations in different languages

C++

/* fft.cpp

* 
* This is a KISS implementation of
* the Cooley-Tukey recursive FFT algorithm.
* This works, and is visibly clear about what is happening where.
*
* To compile this with the GNU/GCC compiler:
* g++ -o fft fft.cpp -lm
*
* To run the compiled version from a *nix command line:
* ./fft
*
*/
  1. include <complex>
  2. include <cstdio>
  1. define M_PI 3.14159265358979323846 // Pi constant with double precision

using namespace std;

// separate even/odd elements to lower/upper halves of array respectively. // Due to Butterfly combinations, this turns out to be the simplest way // to get the job done without clobbering the wrong elements. void separate (complex<double>* a, int n) {

   complex<double>* b = new complex<double>[n/2];  // get temp heap storage
   for(int i=0; i<n/2; i++)    // copy all odd elements to heap storage
       b[i] = a[i*2+1];
   for(int i=0; i<n/2; i++)    // copy all even elements to lower-half of a[]
       a[i] = a[i*2];
   for(int i=0; i<n/2; i++)    // copy all odd (from heap) to upper-half of a[]
       a[i+n/2] = b[i];
   delete[] b;                 // delete heap storage

}

// N must be a power-of-2, or bad things will happen. // Currently no check for this condition. // // N input samples in X[] are FFT'd and results left in X[]. // Because of Nyquist theorem, N samples means // only first N/2 FFT results in X[] are the answer. // (upper half of X[] is a reflection with no new information). void fft2 (complex<double>* X, int N) {

   if(N < 2) {
       // bottom of recursion.
       // Do nothing here, because already X[0] = x[0]
   } else {
       separate(X,N);      // all evens to lower half, all odds to upper half
       fft2(X,     N/2);   // recurse even items
       fft2(X+N/2, N/2);   // recurse odd  items
       // combine results of two half recursions
       for(int k=0; k<N/2; k++) {
           complex<double> e = X[k    ];   // even
           complex<double> o = X[k+N/2];   // odd
                        // w is the "twiddle-factor"
           complex<double> w = exp( complex<double>(0,-2.*M_PI*k/N) );
           X[k    ] = e + w * o;
           X[k+N/2] = e - w * o;
       }
   }

}

// simple test program int main () {

   const int nSamples = 64;
   double nSeconds = 1.0;                      // total time for sampling
   double sampleRate = nSamples / nSeconds;    // n Hz = n / second 
   double freqResolution = sampleRate / nSamples; // freq step in FFT result
   complex<double> x[nSamples];                // storage for sample data
   complex<double> X[nSamples];                // storage for FFT answer
   const int nFreqs = 5;
   double freq[nFreqs] = { 2, 5, 11, 17, 29 }; // known freqs for testing
   
   // generate samples for testing
   for(int i=0; i<nSamples; i++) {
       x[i] = complex<double>(0.,0.);
       // sum several known sinusoids into x[]
       for(int j=0; j<nFreqs; j++)
           x[i] += sin( 2*M_PI*freq[j]*i/nSamples );
       X[i] = x[i];        // copy into X[] for FFT work & result
   }
   // compute fft for this data
   fft2(X,nSamples);
   
   printf("  n\tx[]\tX[]\tf\n");       // header line
   // loop to print values
   for(int i=0; i<nSamples; i++) {
       printf("% 3d\t%+.3f\t%+.3f\t%g\n",
           i, x[i].real(), abs(X[i]), i*freqResolution );
   }

}


C#

using System; using System.Numerics;

class FT {

   public Complex[] x;                // storage for sample data
   public Complex[] X;                // storage for FFT answer
    // separate even/odd elements to lower/upper halves of array respectively.
   // Due to Butterfly combinations, this turns out to be the simplest way 
   // to get the job done without clobbering the wrong elements.
   void Separate(ref Complex[] a, int m, int n)
   {
       Complex[] b = new Complex[(n - m) / 2];
       for (int i = 0; i < (n - m) / 2; i++)    // copy all odd elements to b
           b[i] = a[m + i * 2 + 1];
       for (int i = 0; i < (n - m) / 2; i++)    // copy all even elements to lower-half of a
           a[m + i] = a[m + i * 2];
       for (int i = 0; i < (n - m) / 2; i++)    // copy all odd (from b) to upper-half of a[]
           a[m + i + (n - m) / 2] = b[i];
   }
   void FFT2(ref Complex[] X, int m, int n) // m: the first element of the array X
                                            // n - 1: the last element of the array X
   {
       if (n - m < 2)
       {
           // bottom of recursion.
           // Do nothing here, because already X[0] = x[0]
       }
       else
       {
           Separate(ref X, m, n);      // all evens to lower half, all odds to upper half
           FFT2(ref X, m, m + (n - m) / 2);   // recurse even items
           FFT2(ref X, m + (n - m) / 2, n);   // recurse odd  items                                               
           for (int k = 0; k < (n - m) / 2; k++)// combine results of two half recursions
           {
               Complex e = X[m + k];   // even
               Complex o = X[m + k + (n - m) / 2];   // odd                                                      
               Complex w = Complex.Exp(new Complex(0, -2 * Math.PI * k / (n - m))); // w is the "twiddle-factor"
               X[m + k] = e + w * o;
               X[m + k + (n - m) / 2] = e - w * o;
           }
       }
   }
   double signal(double t)
   {
       double[] freq = { 2, 5, 11, 17, 29 }; // known freqs for testing
       double sum = 0;
       for (int j = 0; j < freq.GetLength(0); j++) // sum several known sinusoids into x[]
           sum += Math.Sin(2 * Math.PI * freq[j] * t);
       return sum;
   }
   void Fill()
   {
       // generate samples for testing
       int N = x.Length;
       for (int i = 0; i < N; i++)
       {
           x[i] = signal((double)i / N);
           X[i] = x[i];  // copy into X[] for FFT work & result
       }
   }
   public FT(int N)
   {
       x = new Complex[N];                // storage for sample data
       X = new Complex[N];                // storage for FFT answer
       Fill();        
       FFT2(ref X, 0, N);
   }

}

class Program {

   static void Main()
   {
       int N = 64;
       FT ft = new FT(N);
       Console.WriteLine("  n\tx[]\tX[]\tf");
       for (int i = 0; i < N; i++)
       {
           Console.WriteLine(i.ToString() + "\t" + ft.x[i].Real.ToString("0.###") + "\t"
               + Complex.Abs(ft.X[i]).ToString("0.###") + "\t" + i);
       }
       Console.ReadKey();
   }

}

An Example

https://github.com/kosme/arduinoFFT/tree/master/Examples

Authors

  • Jordan Gewirtz
  • Nish Chakraburtty
  • Chanel Lynn

Group Link

eBox project page here

eBox log page here

External References

https://www.youtube.com/watch?v=spUNpyF58BY&t=21s

https://www.algorithm-archive.org/contents/cooley_tukey/cooley_tukey.html

https://www.youtube.com/watch?v=XtypWS8HZco&t=22s

http://people.scs.carleton.ca/~maheshwa/courses/5703COMP/16Fall/FFT_Report.pdf

https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm