Go Down

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

deif

Oct 10, 2010, 03:42 pm Last 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[3]);
       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[256]);

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 pm Last 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[128];
char data[128];

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 pm Last 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!
When the testing is complete there will be... cake.

Rusty in Texas

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

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
wr =  pgm_read_word_near(Sinewave + j+N_WAVE/4);


Thanks.

Go Up