Go Down

### Topic: Modified 8bit FFT in c (Read 146350 times)previous topic - next topic

#### deif ##### Oct 10, 2010, 03:42 pmLast Edit: Oct 10, 2010, 03:44 pm by deif Reason: 1
I was looking for a simple arduino fft library to create an equalizer.
I found lots of assembler code or 16 or 32 bit versions but nothing that was simple and out of the box usable.
I have modified the fix_fft.c to use 8bit values.
here is it for everyone with the same problem.

I did not test the inverse fft. I guess there is a bug with the scaling.

fix_fft.h
Code: [Select]
`#ifndef FIXFFT_H#define FIXFFT_H#include <WProgram.h>/*  fix_fft() - perform forward/inverse fast Fourier transform.  fr[n],fi[n] are real and imaginary arrays, both INPUT AND  RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to  0 for forward transform (FFT), or 1 for iFFT.*/int fix_fft(char fr[], char fi[], int m, int inverse);/*  fix_fftr() - forward/inverse FFT on array of real numbers.  Real FFT/iFFT using half-size complex FFT by distributing  even/odd samples into real/imaginary arrays respectively.  In order to save data space (i.e. to avoid two arrays, one  for real, one for imaginary samples), we proceed in the  following two steps: a) samples are rearranged in the real  array so that all even samples are in places 0-(N/2-1) and  all imaginary samples in places (N/2)-(N-1), and b) fix_fft  is called with fr and fi pointing to index 0 and index N/2  respectively in the original array. The above guarantees  that fix_fft "sees" consecutive real samples as alternating  real and imaginary samples in the complex array.*/int fix_fftr(char f[], int m, int inverse);#endif`

fix_fft.cpp
Code: [Select]
`#include <avr/pgmspace.h>#include "fix_fft.h"#include <WProgram.h>/* fix_fft.c - Fixed-point in-place Fast Fourier Transform  *//*  All data are fixed-point short integers, in which -32768  to +32768 represent -1.0 to +1.0 respectively. Integer  arithmetic is used for speed, instead of the more natural  floating-point.  For the forward FFT (time -> freq), fixed scaling is  performed to prevent arithmetic overflow, and to map a 0dB  sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq  coefficients. The return value is always 0.  For the inverse FFT (freq -> time), fixed scaling cannot be  done, as two 0dB coefficients would sum to a peak amplitude  of 64K, overflowing the 32k range of the fixed-point integers.  Thus, the fix_fft() routine performs variable scaling, and  returns a value which is the number of bits LEFT by which  the output must be shifted to get the actual amplitude  (i.e. if fix_fft() returns 3, each value of fr[] and fi[]  must be multiplied by 8 (2**3) for proper scaling.  Clearly, this cannot be done within fixed-point short  integers. In practice, if the result is to be used as a  filter, the scale_shift can usually be ignored, as the  result will be approximately correctly normalized as is.  Written by:  Tom Roberts  11/8/89  Made portable:  Malcolm Slaney 12/15/94 malcolm@interval.com  Enhanced:  Dimitrios P. Bouras  14 Jun 2006 dbouras@ieee.org  Modified for 8bit values David Keller  10.10.2010*/#define N_WAVE      256    /* full length of Sinewave[] */#define LOG2_N_WAVE 8      /* log2(N_WAVE) *//*  Since we only use 3/4 of N_WAVE, we define only  this many samples, in order to conserve data space.*/const prog_int8_t Sinewave[N_WAVE-N_WAVE/4] PROGMEM = {0, 3, 6, 9, 12, 15, 18, 21, 24, 28, 31, 34, 37, 40, 43, 46, 48, 51, 54, 57, 60, 63, 65, 68, 71, 73, 76, 78, 81, 83, 85, 88, 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 109, 111, 112, 114, 115, 117, 118, 119, 120, 121, 122, 123, 124, 124, 125, 126, 126, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 127, 126, 126, 125, 124, 124, 123, 122, 121, 120, 119, 118, 117, 115, 114, 112, 111, 109, 108, 106, 104, 102, 100, 98, 96, 94, 92, 90, 88, 85, 83, 81, 78, 76, 73, 71, 68, 65, 63, 60, 57, 54, 51, 48, 46, 43, 40, 37, 34, 31, 28, 24, 21, 18, 15, 12, 9, 6, 3, 0, -3, -6, -9, -12, -15, -18, -21, -24, -28, -31, -34, -37, -40, -43, -46, -48, -51, -54, -57, -60, -63, -65, -68, -71, -73, -76, -78, -81, -83, -85, -88, -90, -92, -94, -96, -98, -100, -102, -104, -106, -108, -109, -111, -112, -114, -115, -117, -118, -119, -120, -121, -122, -123, -124, -124, -125, -126, -126, -127, -127, -127, -127, -127, /*-127, -127, -127, -127, -127, -127, -126, -126, -125, -124, -124, -123, -122, -121, -120, -119, -118, -117, -115, -114, -112, -111, -109, -108, -106, -104, -102, -100, -98, -96, -94, -92, -90, -88, -85, -83, -81, -78, -76, -73, -71, -68, -65, -63, -60, -57, -54, -51, -48, -46, -43, -40, -37, -34, -31, -28, -24, -21, -18, -15, -12, -9, -6, -3, */};/*  FIX_MPY() - fixed-point multiplication & scaling.  Substitute inline assembly for hardware-specific  optimization suited to a particluar DSP processor.  Scaling ensures that result remains 16-bit.*/inline char FIX_MPY(char a, char b){    //Serial.println(a); //Serial.println(b);        /* shift right one less bit (i.e. 15-1) */    int c = ((int)a * (int)b) >> 6;    /* last bit shifted out = rounding-bit */    b = c & 0x01;    /* last shift + rounding bit */    a = (c >> 1) + b;        /*        Serial.println(Sinewave);        Serial.println(c);        Serial.println(a);        while(1);*/    return a;}/*  fix_fft() - perform forward/inverse fast Fourier transform.  fr[n],fi[n] are real and imaginary arrays, both INPUT AND  RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to  0 for forward transform (FFT), or 1 for iFFT.*/int fix_fft(char fr[], char fi[], int m, int inverse){    int mr, nn, i, j, l, k, istep, n, scale, shift;    char qr, qi, tr, ti, wr, wi;    n = 1 << m;    /* max FFT size = N_WAVE */    if (n > N_WAVE)        return -1;    mr = 0;    nn = n - 1;    scale = 0;    /* decimation in time - re-order data */    for (m=1; m<=nn; ++m) {        l = n;        do {            l >>= 1;        } while (mr+l > nn);        mr = (mr & (l-1)) + l;        if (mr <= m)            continue;        tr = fr[m];        fr[m] = fr[mr];        fr[mr] = tr;        ti = fi[m];        fi[m] = fi[mr];        fi[mr] = ti;    }    l = 1;    k = LOG2_N_WAVE-1;    while (l < n) {        if (inverse) {            /* variable scaling, depending upon data */            shift = 0;            for (i=0; i<n; ++i) {                j = fr[i];                if (j < 0)                    j = -j;                m = fi[i];                if (m < 0)                    m = -m;                if (j > 16383 || m > 16383) {                    shift = 1;                    break;                }            }            if (shift)                ++scale;        } else {            /*              fixed scaling, for proper normalization --              there will be log2(n) passes, so this results              in an overall factor of 1/n, distributed to              maximize arithmetic accuracy.            */            shift = 1;        }        /*          it may not be obvious, but the shift will be          performed on each data point exactly once,          during this pass.        */        istep = l << 1;        for (m=0; m<l; ++m) {            j = m << k;            /* 0 <= j < N_WAVE/2 */            wr =  pgm_read_word_near(Sinewave + j+N_WAVE/4);/*Serial.println("asdfasdf");Serial.println(wr);Serial.println(j+N_WAVE/4);Serial.println(Sinewave);Serial.println("");*/            wi = -pgm_read_word_near(Sinewave + j);            if (inverse)                wi = -wi;            if (shift) {                wr >>= 1;                wi >>= 1;            }            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;            }        }        --k;        l = istep;    }    return scale;}/*  fix_fftr() - forward/inverse FFT on array of real numbers.  Real FFT/iFFT using half-size complex FFT by distributing  even/odd samples into real/imaginary arrays respectively.  In order to save data space (i.e. to avoid two arrays, one  for real, one for imaginary samples), we proceed in the  following two steps: a) samples are rearranged in the real  array so that all even samples are in places 0-(N/2-1) and  all imaginary samples in places (N/2)-(N-1), and b) fix_fft  is called with fr and fi pointing to index 0 and index N/2  respectively in the original array. The above guarantees  that fix_fft "sees" consecutive real samples as alternating  real and imaginary samples in the complex array.*/int fix_fftr(char f[], int m, int inverse){    int i, N = 1<<(m-1), scale = 0;    char tt, *fr=f, *fi=&f[N];    if (inverse)        scale = fix_fft(fi, fr, m-1, inverse);    for (i=1; i<N; i+=2) {        tt = f[N+i-1];        f[N+i-1] = f[i];        f[i] = tt;    }    if (! inverse)        scale = fix_fft(fi, fr, m-1, inverse);    return scale;}`

#### deif #1
##### Oct 10, 2010, 03:43 pmLast Edit: Oct 10, 2010, 04:05 pm by deif Reason: 1
this is the way I used it:
have fun sketch:
Code: [Select]
`char im;char data; void loop(){  int static i = 0;  static long tt;  int val;     if (millis() > tt){      if (i < 128){        val = analogRead(pin_adc);        data[i] = val / 4 - 128;        im[i] = 0;        i++;                }      else{        //this could be done with the fix_fftr function without the im array.        fix_fft(data,im,7,0);        // I am only interessted in the absolute value of the transformation        for (i=0; i< 64;i++){           data[i] = sqrt(data[i] * data[i] + im[i] * im[i]);        }                //do something with the data values 1..64 and ignore im        show_big_bars(data,0);      }        tt = millis();   }}`

this should give an fft with
- sampling rate:           1ms
- frequency resolution:  500Hz
- lowest frequency:       7.8Hz

there is a lot that can be improved:
- fix bugs
- replace FIX_MPY by inline assembler
- do sampling with timer and interrupts
- use fix_fftr instead of fix_fft to save ram (i did not get the math right)

#### focalist #2
##### Nov 14, 2010, 06:30 pmLast Edit: Nov 14, 2010, 06:54 pm by focalist Reason: 1
Hmm.   Okay, here's the good news:  I've used your code, and it works wonderfully for my purposes.  Will be posting a project doc soon on it.  This code is probably going to be useful for a number of projects, much praise for doing it, and even more for putting it up here for general use.

The "bad" news is that I'm still a little confused about the function call parameter "m".  Note that I am not a mathematics expert- I've simply used your code as a black box audio-range FFT with 64 band resolution, just as you have.  I then diddle with outputting the array in various ways.  Works great, just working out gain on amplifier for electret mic.  It works with a single-stage NPN, but the signals are pretty small.. later today going to add a second NPN (collector direct-coupled, E1-B2 coupled, with and without capacitor, to make a darlington of them) to see if the signal is better.. would be a nice thing if works out to be that easy (doubtful).

I'm trying to understand the underlying math, it's some pretty hardcore stuff.. matrix math and imaginary components.. erm... where's the aspirin? So instead, I'm starting with examining what you have done, to try to get a mental handle on the process from a practical standpoint- as the implementation works, I want to understand WHY it works.

I've got my head around the arrays and why they are needed, but what exactly is "m", in practical terms?  It's cryptically referred to as a bit-shift parameter for the data, but that seems wrong too.  Is it number of data bits (7 bits, 0-127) in the output array data?

Thanks for any insight, I'm amazed our little arduino can perform a transform reasonably... this code is perfect for my uses, so I want to understand at least the function call parameters, even if the rest remains a Black Box for the time being..

Great work!

#### Fe2o3Fish #3
##### Nov 15, 2010, 07:32 am
m = log2(n) where n is the number of data points in the input.

log2(n) =  the logarithm, base 2, of n. Essentially the number of
bits in the value of n.

Since the size of data[] is 128, the n is 7 because 2**7 == 128.

The number of input points/samples for FFT's should be a nice
power of 2 (8,16,32,etc).  If not, zero-filled at the end to fill that
data array to the next highest power of two before an FFT is
performed.

Hope that helps.

#### tendor #4
##### Nov 16, 2010, 08:17 am
Hello, deif.
You use int8_t ie 8bit char here
Quote
const prog_int8_t Sinewave[N_WAVE-N_WAVE/4] PROGMEM = {

And why do you read words into the fix_fft() function?
Quote

Thanks.

#### Cabe #5
##### Nov 19, 2010, 10:46 pm
Musical christmas light and DIY DJ's across the globe now have a new friend.

Also featured on Make - http://blog.makezine.com/archive/2010/11/real-time_spectrum_analyzer_powered.html

#### sydkahn #6
##### Nov 26, 2010, 11:20 am
Calculus not being my fote - can the FFT be broken/summed up into 4 frequency blocks?

I have a few of http://diylightanimation.com/wiki/index.php?title=Manual_SSR4- that are four channel SSR type mains light controllers - and would like break down into 4 channels to drive each switch.

I built scads of Color organs back in the late 60's, so i understand how to turn on a triac, but doing it digitally seems to be lots more fun...

#### focalist #7
##### Nov 27, 2010, 11:12 amLast Edit: Nov 27, 2010, 11:15 am by focalist Reason: 1
Thanks, Rusty... been reading the code and thought that was the case, thanks for the confirmation.

Syd, that would be very easy... what you are going to do is average bands together.  The build I did with TvOut uses 64 bands.. I just display all 64 bands (0-32kHz approximately, 500Hz/band).  What you want is to reduce the number of bands, and that's done by averaging across the data.  Since you want four bands out, that means you want to take the first 16 data points, add 'em together, and divide by 16.  That's your value for band #1, which would be 0-8Khz.  The next 16 data points are the second band, etc.

This will give you four values, to do with as you please.  Since you are using simple on/off, do a bit of testing (or put four potentiometers on the Arduino's A1-A4) to set a threshhold value to switch the SSR's on (light up when sound in that band is above threshhold level), send that on/off (High/Low, 1/0) to four digital pins, hook up your SSR's to those pins.  If you have them around, the Pots are probably the easiest way to set threshholds as real-world input levels (and desired response) change.

Working on an updated version of this, now that Turkey Day is past..

#### focalist #8
##### Nov 27, 2010, 03:14 pmLast Edit: Nov 27, 2010, 03:16 pm by focalist Reason: 1
Well Syd, you're going to be as shocked as I was with how fast the FFT runs.. in my implementation, not only is our little Arduino doing the sampling and FFT's, it's also generating composite video, completely in code, thanks to the TVout library.  I'm still doing some tweaking (after which I will do actual timing to see the exact rate), but with the build I wrote up, I'm guessing it's updating (which means a sample/process/display pass) the screen at around 10 times per second.  A 10Hz frame rate isn't going to be too good for an analytical device, but for a light show, it'll do just fine.  The majority of the time is actually being consumed by the line-drawing routines for the bars, rather than the data transform math... since your code handling output will no doubt be less "weighty" than generating NTSC video, you will probably get a MUCH higher update rate.

#### focalist #9
##### Nov 28, 2010, 06:45 pmLast Edit: Nov 28, 2010, 06:59 pm by focalist Reason: 1
By the way, Arduino has done it to me once again... a project is forcing me to learn more.  My successful use of this code made me want to understand it better- wasn't happy with a black box, I like to see the inner workings and understand them if I can.

So, I've been reading through the code.. and then wiki.. and then EE society writeups.. and university papers... etc etc etc.  The whole conversion to integer math makes the understanding even more difficult.. bit shift math aside, it's way to easy to get confused with the data format.

In any case, I'm now learning about the math behind this, because of you... dangit... lol.. Arduino strikes again!

The best part is that it is OPEN SOURCE.  This allowed me to use the original work without needing to understand the complexities behind the scenes- but when I decided I wanted to see the inner workings, they are right there for the reading.. often with well-commented code and external references by the author for THEIR foundations.  Simply put, you can pick it apart to the lowest level if you desire.  Same with the TVout lib.  These open source libs, on an open source platform- they allowed me to build a project far beyond my own technical abilities with relative ease.  This kind of stuff is what Open Source and the Arduino project are really all about, in my humble opinion.  If given enough time and reference material, could I have produced the project on my own?  Maybe.  Doubtful, in fact pretty unlikely.  Even if I did, the effort level involved would just be too much to derive "fun" from it.  Months of work, hardcore stuff.  With Arduino (and the community), I went from idea to working prototype in a couple of hours.  How can you beat that?  I'll tell you how:  by putting the source code out, it allows me to BACKTRACK to understand the parts I want to.. AFTER I've shown myself the idea works...

#### DE8MSH #10
##### Nov 30, 2010, 09:56 pmLast Edit: Nov 30, 2010, 10:16 pm by DE8MSH Reason: 1
Hi,

hope that anyone can help: will use this FFT for some HAM projects. Question: how can I reduce the bandwidth from 330hz to 3300 khz only?

Most HAM radios have audio bandwidths from 0 to 3khz.

And: how can I use 128 or 256 FFT level instead of 64?

Thank you very much fpr information.

#### RicardoLATech #11
##### Dec 06, 2010, 01:08 am
I'm currently trying to take an analog signal into my arduino in the frequency range of around 70-400 Hz with a few harmonics in each signal.  Can this FFT code be used if i wanted to find the fundamental frequency of my signals?  How would i go about doing that?  Do i need to look for the highest magnitude in the FFT bins?  I'm new to this stuff and need a little guidance.

Thanks,

#### focalist #12
##### Dec 07, 2010, 04:27 pm
Ricardo, yes.. that's going to be your best bet.  There's a bit of error in the bands because of the use of Integer math, but it's at the cost of speed to go to floating point of any sort.

Unlike my above project, you are going to want to actually tune the sample code to actually get timed samples, as my approximation (or even "fast as you can" code) is NOT lab accurate by any means.

I've been reading up and am starting to get a grip on the math.  By the way, to the other poster.. the output is N/2 where N is the number of samples, this gives 64 bins due to 128 byte sample space.

I might tweak this code at some point and stick it somewhere as a downloadable lib-- I haven't seen the author back to the thread, but this code shouldn't disappear.. far too useful, if even just for a light show..

#### NoX #13
##### Jan 04, 2011, 07:46 pm
Hi all!

I have been trying to get this working with no luck, I have tried versions 18, 21 and 22 in ubuntu 64bits and in windows xp

I allways get this error:
Code: [Select]
`FFT_test.cpp.o: In function `setup':FFT_test.cpp:21: undefined reference to `fix_fft(char*, char*, int, int)'`

I got the code for the pde from http://blurtime.blogspot.com/2010/11/arduino-realtime-audio-spectrum.html
but I think the problem is in the library... 