Fast Fourier Transform Library?

Is there an FFT library available? I’d like to hook up a microphone and analyze sound data.

3 Likes

http://www.arduinoos.com/2010/10/fast-fourier-transform-fft/

The Teensy library looks like a good starting point considering it’s using a 32-bit ARM Cortex-M4 processor running at 48MHz

1 Like

dougal, BDub’s reference to the CMSIS library used in the teensy3 code is fantastic. There is also this arduino library you could portt:

http://wiki.openmusiclabs.com/wiki/ArduinoFFT

If you want to punish yourself, here is some Cortex M3 assembler libraries:
http://www.embeddedsignals.com/ARM.htm

:smile:

2 Likes

Given that space on the Spark is in short supply but the processor is quite brisk I found that this very simple code worked …

/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
#include <math.h>

short FFT(short int dir,int m,float *x,float *y)
{
   int n,i,i1,j,k,i2,l,l1,l2;
   float c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }
   
   return(0);
}

I think that the open music labs code uses low level Atmega code that won’t port easily.

3 Likes

Thanks for all the links, everybody. I skimmed through the Arduinoos.com links to refresh my memory on how to actually put FFT to use. I’ll try playing with it when I have time.

Or who knows, for what I want to do, I might decide to just implement some bandpass filters in hardware and read the results on analog pins.

@dougal, from your opening post I gathered that you’d use FFT to analyze a mic input.
FFT for one thing should not be the big issue, given the whizzy processor (as prev. posts show), but how have you planned to hook the mic to the Core?
As far as I understood there are several problems to overcome if you want to use the onboard ADCs (e.g. Cloud heartbeat, ADC speed, sampling frequency, …)

Plan? What is this plan you speak of? :wink:

I'm just going to experiment. I have an Adafruit electret mic w/pre-amp that I'm going to play with. But as for the rest, I'm just winging it.

2 Likes

I thought this sounded like a plan :wink:


:thought_balloon: Planning (also called forethought) is the process of thinking about and organizing the activities required to achieve a desired ("I'd like to ...") goal

@dougal, Apart from CMSIS provided DSP libs, ST also has a library for STM32 implemented in ASM. I have created a new features branch in core-common-lib with the new DSP library additions.

For IIR filter design, my personal preference is mkfilter design for coefficient generation if MATLAB is not available: http://www-users.cs.york.ac.uk/~fisher/mkfilter/

The header file “stm32-dsp.h” contains all the needed functions.

5 Likes

That is quite nice as a free tool. Things to look out for:

  • If you ask for a high-order IIR, you are almost always better off stability and numeric-wise with a second-order section filter where the larger polynomial is broken up into many second-order polynomials. MATLAB can do some optimal reordering of the sections which can really help out in some cases, particularly with Chebyshev filters.
  • The generated C code evaluates the filter polynomial directly but you would be better off using Horner's method in floating point.

I wanted to also say that the nice in-place FFT presented above does scaling in a way that differs from MATLAB. With fft and ifft transforms, there is a gain of N if you just run them back to back, so to make a transform pair, you must have a gain of 1/N somewhere. In MATLAB, the 1/N goes on the ifft. There is no right or wrong answer to where to place this gain--it is just convention, but most signal processing textbooks put the 1/N gain in the inverse side. Other fields like thermodynamics and economics put it on the forward transform.

There are a bunch of tricks if you are doing a real-only input FFT. The Rick Lyons book is a nice reference.

1 Like

Now we can also use the kiss_fft routine from this 3D cube project
kiss fft

Brahim

@Phec – would it be possible for you to post your complete code that uses this FFT function? I’m tinkering with some audio code, and I could use a complete example to help me understand how to get started with your code

Thanks!
–alex

@enjrolas try something like that ...

Note 1. the results looks similar to MATLAB
Note 2. @Phec routine looks like it is from Fast Fourier Transform
See Appendix B. FFT (Fast Fourier Transform)

short int dir;
int m;
float x[8]; // 2^3 = 8
float y[8]; // 2^3 = 8

dir = 1;
m = 3;

x[0] = 0.0;
x[1] = 0.1;
x[2] = 0.2;
x[3] = 0.3;
x[4] = 0.4;
x[5] = 0.3;
x[6] = 0.2;
x[7] = 0.1;

y[0] = 0.0;
y[1] = 0.0;
y[2] = 0.0;
y[3] = 0.0;
y[4] = 0.0;
y[5] = 0.0;
y[6] = 0.0;
y[7] = 0.0;

FFT(dir, m, (float*)x, (float*)y);

// One cane use some printf to check and debug the results 
printf("x[0] = %10.15f \n", x[0]);
printf("x[1] = %10.15f \n", x[1]);
... 
printf("y[0] = %10.15f \n", y[0]);
printf("y[1] = %10.15f \n", y[1]);

@bhamadicharef Thank you for the reference to the source of the FFT - my code was a copy of a copy so I’m very pleased that the original author is now identified.
Alex, @enjrolas - the code I posted works in a short stand alone program on the Spark or Arduino (where it is pinched from) however there isn’t a lot of RAM to play with so I haven’t used the function on the Spark in anger - I do my FFT on a Raspberry Pi. The PI python code is in the [Beehive Monitor Git Gist ][1].

However if you don’t have too big an array so that memory isn’t an issue here is a short program that calls the FFT function from my earlier post. I haven’t used any windowing:

float x[16];
float y[16];
// copy the FFT function here
void setup(){
      //put a sine wave into x 
      for (int i=0;i<16;i++){
        x[i]=sin(1.0*i*3.142/8);
        y[i]=0;       
    }
    
    Serial.begin(9600);
    for (int i=0;i<16;i++){
      Serial.println(sqrt(x[i]*x[i]+y[i]*y[i]));
    }
// do the FFT 1=forward, 4=16 bins, x and y contain the input and output
    FFT(1,4,x,y);
    for (int i=0;i<16;i++){
      Serial.println(y[i]);
    }
}

void loop() {

}






I have a couple of photons on order and since they have more generous amounts of RAM I plan to put an fft on them so I’m interested in the development of an efficient Spark fft library.
[1]: https://gist.github.com/phec/23cb3cfb27c955e182e6#file-beemonitorblackufl14-py

1 Like

@bhamadicharef -- Thanks! It made sense to me after reading the paul bourke page. I hadn't understood the "transform in place" idea. Awesome!

@satishgn Thanks for posting that DSP library! How can I include it in my code if it is not in the main branch?

So I was reviewing my old posts, and this piqued my interest again… I ran across a post about a SplitRadix library. (in particular, the SplitRadixRealP code, for the Arduino DUE)

This library is supposed to perform FFT with much less RAM usage and greater speed than most other implementations. The library itself seems to compile okay in the web IDE. But the demo code does not, as it contains lots of PIO_* and TC_* references which aren’t recognized.

I haven’t tried to decipher what that code is doing, myself. Does anyone want to take a stab at porting it (@peekay123?) ? Or is it even necessary to do so?

Hello together,

i am trying to achieve something similar. I have an issue with the arm math library. Maybe you can help me to move forward.

https://community.particle.io/t/arm-math-perform-fft/44314

Thank you!