8-bit FFT revisited (FIX_FFT)

It rained (again) in Georgia, so my wife left me alone in my lab to play yesterday XD
I have been playing around with an Atmega1284P-PU on one of CrossRoad's boards; for some reason, I decided to try the FFT algorithm. Searching the Forum and Google was interesting (not in an overly positive way), so when I finally got something to work based on an old forum post: Modified 8bit FFT in c - Development - Arduino Forum I decided to clean-up the code a bit, add some small comments, and double-check the performance with my Rigol oscilloscope. The results of that is in this post.

The last entry by Magician back in January of '11 was

Unfortunately, there is a flow in the code, no FFT performed
if you copy/paste as it is.
:-[

My goal for this post was to fix this cut & paste issue.

I have checked the following test code on UNO and Bobuino_1284 and it works fine using Arduino 1.0.5. To avoid having to set up the .h and the .cpp files as libraries, just put them into your Arduino sketch folder with the .INO file. I do this to allow me to make minor changes (ONLY IF Necessary) before migrating same to the standard \Arduino\Libraries folder. It is just a preference of mine and certainly is not necessary, but if you use the traditional approach, be sure to change #include "fix_fft.h" to #include <fix_fft.h> so that the GUI can find the files at compile time.

Based on my testing, (I do not use bin #0, rather using the display space for the "L" and "R" indicating the Left and Right channels) bin#1 is tuned to 85Hz and the following 14 bins are: 120, 190, 240, 305, 365, 430, 480, 550, 625, 660, 735, 780, 850, and 910Hz. You will note that it is not exactly linear because we are not using an interrupt-driven paradigm for accurate timing; that is, the Analog inputs are just read in a loop. Accurate timing examples are given in some of the music note implementations of FFT use.

Audio signal is approximately 1V P-P from the audio generator, fed through a 0.5uF (0.5MFD value not critical) to the junction of two (2) 10K resistors with the other ends of the resistors going to +5V and Gnd. This voltage-divider junction establishes the 2.5V balance to the AD port. Note, you must use a capacitor and two resistors for each of the two analog inputs. In my breadboard design, I used a 680 Ohm resistor between the voltage divider junction and the input to the uC analog port as a kind of safety to limit any mishap to around 7mA of current BUT this is not necessary in the final layout unless you just have a bunch of 680-1000 Ohm resistors laying around (I do.)

Main code:

/* FFT_TEST4
 Ray Burnette 20130810 function clean-up & 1284 port (328 verified)
 Uses 2x16 Parallel LCD in 4-bit mode, see LiquidCrystal lib call for details
 http://forum.arduino.cc/index.php?PHPSESSID=4karr49jlndufvtlqs9pdd4g96&topic=38153.15
 Modified by varind in 2013: this code is public domain, enjoy!
 http://www.variableindustries.com/audio-spectrum-analyzer/
 328P = Binary sketch size: 5,708 bytes (of a 32,256 byte maximum)
 1284P= Binary sketch size: 5,792 bytes (of a 130,048 byte maximum) Free RAM = 15456
        Binary sketch size: 8,088 bytes (of a 130,048 byte maximum) (Debug)
 */

#include <LiquidCrystal.h>
#include "fix_fft.h"  // fix_fft.ccp & fix_fft.h in same directory as sketch

#define DEBUG 0
#define LCHAN 1
#define RCHAN 0

const int Yres = 8;
const int gain = 3;
float peaks[64];
char im[64], data[64];
char Rim[64], Rdata[64];
char data_avgs[64];
int debugLoop;

//LiquidCrystal(rs, enable, d4, d5, d6, d7) 1284P Physical: 6,  5, 4, 3, 2, 1
LiquidCrystal lcd(11, 10, 7, 6, 5, 4); // saves all analog pins port PA

// Custom CHARACTERS
byte v1[8] = {
  B00000,B00000,B00000,B00000,B00000,B00000,B00000,B11111};
byte v2[8] = {
  B00000,B00000,B00000,B00000,B00000,B00000,B11111,B11111};
byte v3[8] = {
  B00000,B00000,B00000,B00000,B00000,B11111,B11111,B11111};
byte v4[8] = {
  B00000,B00000,B00000,B00000,B11111,B11111,B11111,B11111};
byte v5[8] = {
  B00000,B00000,B00000,B11111,B11111,B11111,B11111,B11111};
byte v6[8] = {
  B00000,B00000,B11111,B11111,B11111,B11111,B11111,B11111};
byte v7[8] = {
  B00000,B11111,B11111,B11111,B11111,B11111,B11111,B11111};
byte v8[8] = {
  B11111,B11111,B11111,B11111,B11111,B11111,B11111,B11111};


void setup() {

  if (DEBUG) {
    Serial.begin(9600); // hardware serial
    Serial.print("Debug ON");
    Serial.println("");
  }

  lcd.begin(16, 2);
  lcd.clear();
  lcd.createChar(1, v1);
  lcd.createChar(2, v2);
  lcd.createChar(3, v3);
  lcd.createChar(4, v4);
  lcd.createChar(5, v5);
  lcd.createChar(6, v6);
  lcd.createChar(7, v7);
  lcd.createChar(8, v8);
}

void loop() {

  for (int i = 0; i < 64; i++) {    // 64 bins = 32 bins of usable spectrum data
    data[i]  = ((analogRead(LCHAN) / 4 ) - 128);  // chose how to interpret the data from analog in                                      
    im[i]  = 0;   // imaginary component
    Rdata[i] = ((analogRead(RCHAN) / 4 ) - 128);  // chose how to interpret the data from analog in                                      
    Rim[i] = 0;   // imaginary component
  }

  fix_fft(data, im, 6, 0);   // Send Left channel normalized analog values through fft
  fix_fft(Rdata, Rim, 6, 0); // Send Right channel normalized analog values through fft

  // At this stage, we have two arrays of [0-31] frequency bins deep [32-63] duplicate

  // calculate the absolute values of bins in the array - only want positive values
  for (int i = 0; i < 32; i++) {
    data[i] = sqrt(data[i]  *  data[i] +  im[i] *  im[i]);
    Rdata[i] = sqrt(Rdata[i] * Rdata[i] + Rim[i] * Rim[i]);

    // COPY the Right low-band (0-15) into the Left high-band (16-31) for display ease
    if (i < 16) {
      data_avgs[i] = data[i];
    } 
    else {
      data_avgs[i] = Rdata[i - 16];
    }

    // Remap values to physical display constraints... that is, 8 display custom character indexes + "_"
    data_avgs[i] = constrain(data_avgs[i], 0, 9 - gain);     //data samples * range (0-9) = 9
    data_avgs[i] = map(data_avgs[i], 0, 9 - gain, 0, Yres);  // remap averaged values
  }

  Two16_LCD();
  decay(1);
}


void Two16_LCD(){
  lcd.setCursor(0, 0); 
  lcd.print("L"); // Channel ID replaces bin #0 due to hum & noise
  lcd.setCursor(0, 1); 
  lcd.print("R"); // ditto

  for (int x = 1; x < 16; x++) {  // init 0 to show lowest band overloaded with hum
    int y = x + 16; // second display line
    if (data_avgs[x] > peaks[x]) peaks[x] = data_avgs[x];
    if (data_avgs[y] > peaks[y]) peaks[y] = data_avgs[y];

    lcd.setCursor(x, 0); // draw first (top) row Left
    if (peaks[x] == 0) {
      lcd.print("_");  // less LCD artifacts than " "
    }
    else {
      lcd.write(peaks[x]);
    }

    lcd.setCursor(x, 1); // draw second (bottom) row Right
    if (peaks[y] == 0){
      lcd.print("_");
    }
    else {
      lcd.write(peaks[y]);
    }
  }
  
  debugLoop++;
  if (DEBUG && (debugLoop > 99)) {
    Serial.print( "Free RAM = " );
    Serial.println( freeRam(), DEC);
    Serial.println( millis(), DEC);
    debugLoop = 0;
  } 
}


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


void decay(int decayrate){
  int DecayTest = 1;
  // reduce the values of the last peaks by 1 
  if (DecayTest == decayrate){
    for (int x = 0; x < 32; x++) {
      peaks[x] = peaks[x] - 1;  // subtract 1 from each column peaks
      DecayTest = 0;
    }
  }

  DecayTest++;
}

Continued in next post...


Ray

Continued for completeness, but no changes except the commenting of the old Arduino header file:

//#include <WProgram.h>

from that shown on post:
http://forum.arduino.cc/index.php?topic=38153.0

fix_fft.h

#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

#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;  //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);
		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;
}

A couple of pixs... The LCD is showing a single tone of 780Hz being fed simultaneously to input A0 and A1 on the Atmega1284 (physical chip pins 33 & 34.)

Note:
Forum member Magician has evolved the fix_fft library through several iterations, each producing a better/faster product. Please visit his site for the latest as well as some excellent Arduino projects:
http://fftarduino.blogspot.com/2011_02_01_archive.html

Thanks Magician!

Ray

I am trying this on a teensy 3.2 running at 72 mhz but I am not getting a proper spectrum (I am getting a random spectrum). Please help.