chipKIT® Development Platform

Inspired by Arduino™

Audio Spectrum Analyser FFT Fast Fourier Transform QVGA LCD

Created Sun, 30 Oct 2011 21:58:00 +0000 by serveurperso

serveurperso

Sun, 30 Oct 2011 21:58:00 +0000

Hi,

Another toy for the ChipKit and the 20$ TFT QVGA LCD : It use fix_fft.c functions

[url]http://www.youtube.com/watch?v=yoquqEL8OCc[/url]

#include <CK_ITDB02_Graph16.h>

#define TFTX 320
#define TFTY 240

#define BACKLIGHTPIN 2

#define LOG2_N_WAVE 11                 // 11 bits
#define N_WAVE (1 << LOG2_N_WAVE)      // 2048 dots
#define N_WAVE34 (N_WAVE - N_WAVE / 4) // Sinus lookup table size (3/4 wave)

#define GAIN 32                        // Preamplifier
#define DOWNSPEED 2                    // Fall down rate (pixels per frame)
#define OFFSETX 1                      // 0 to N_WAVE - TFTX

#define OVERSAMPLING 2                 // 0 for 1x
                                       // 1 for 2x
                                       // 2 for 4x

#define RGBSPEED (2 * PI / 50)         // 50 frames for one RGB cycle

#define PI23 (2 * PI / 3)
#define PI43 (4 * PI / 3)

extern uint8_t SmallFont[];

ITDB02 tft(82, 83, 84, 85, ITDB32S);

int16_t Sinewave[N_WAVE34];

int16_t samples[N_WAVE];
volatile bool samplesReady = false;
volatile uint16_t sampleCounter = 0;

void setup() {
 tft.InitLCD(LANDSCAPE);
 tft.setFont(SmallFont);
 tft.clrScr();

 pinMode(BACKLIGHTPIN, OUTPUT);

 digitalWrite(BACKLIGHTPIN, HIGH);

 Serial.begin(115200);

 // ADC setup
 AD1PCFG = ~1;          // PORTB digital, RB0 analog
 AD1CHSbits.CH0SA = 0;  // Channel 0 positive input is AN0
 //AD1CON2bits.VCFG = 1;  // External VREF+ pin
 AD1CON1bits.SSRC = 7;  // Internal counter ends sampling and starts conversion (auto-convert)
 AD1CON1bits.ASAM = 1;  // Sampling begins when SAMP bit is set

 // ADC Conversion Clock Select
 // TPB = (1 / 80MHz) * 2 = 25ns
 // TAD = TPB * (ADCS + 1) * 2 = 25ns * 20 = 500ns
 AD1CON3bits.ADCS = 19;

 // Auto-Sample Time
 // SAMC * TAD = 30 * 500ns = 15uS
 AD1CON3bits.SAMC = 30;

 IPC6bits.AD1IP = 6;    // Analog-to-Digital 1 Interrupt Priority
 IPC6bits.AD1IS = 3;    // Analog-to-Digital 1 Subpriority
 IFS1CLR = 2;           // Clear ADC interrupt flag
 IEC1bits.AD1IE = 1;    // ADC interrupt enable
 AD1CON1bits.ON = 1;    // ADC enable

 for(uint16_t i = 0; i < N_WAVE34; i++)
  Sinewave[i] = int(sin(i * 2 * PI / N_WAVE) * 32767);
}

void loop() {
 static uint16_t i;
 static uint8_t j;
 static int16_t real[N_WAVE];
 static int16_t imag[N_WAVE];
 static uint16_t amplitude[N_WAVE / 2];

 while(!samplesReady);
 samplesReady = false;

 for(i = 0; i < N_WAVE >> OVERSAMPLING; i++) {
  real[i] = samples[i << OVERSAMPLING];
  for(j = 1; j < 1 << OVERSAMPLING; j++) {
   real[i] += samples[(i << OVERSAMPLING) + j];
  }
  real[i] = (real[i] >> OVERSAMPLING) * GAIN;
  imag[i] = 0;
 }
/*
 for(i = 0; i < N_WAVE >> OVERSAMPLING; i++) {
  real[i] = samples[i << OVERSAMPLING] * GAIN;
  imag[i] = 0;
 }
*/

 AD1CON1bits.ASAM = 1; // Start autosampling

 fix_fft(real, imag, LOG2_N_WAVE, 0);

 for(i = OFFSETX; i < OFFSETX + TFTX; i++)
  amplitude[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);

 drawSpectrum(amplitude);
}

void drawSpectrum(uint16_t amplitude[]) {
 static uint16_t oldY[TFTX] = {0};
 static float t = 0;
 uint16_t y;
 uint16_t yDown;
 uint8_t r = int(cos(t) * 127 + 127);
 uint8_t g = int(cos(t + PI43) * 127 + 127);
 uint8_t b = int(cos(t + PI23) * 127 + 127);

 t += RGBSPEED;
 if(t >= 2 * PI) t = 0;

 for(uint16_t x = 0; x < TFTX; x++) {
  y = amplitude[x + OFFSETX];
  if(y > TFTY) y = TFTY;

  if(y > 0) {
   tft.setColor(r, g, b);
   tft.drawLine(x, TFTY - y, x, TFTY);
  }

  if(y < oldY[x]) {
   if(oldY[x] >= DOWNSPEED)
    yDown = DOWNSPEED;
   else
    yDown = oldY[x];

   tft.setColor(0, 0, 0);
   tft.drawLine(x, TFTY - oldY[x],
                x, TFTY - (oldY[x] - yDown));

   oldY[x] -= yDown;
  } else if(y > 0)
   oldY[x] = y;

 }
}

extern "C" {
 void __ISR(_ADC_VECTOR, ipl6) ADCInterruptHandler() {
  IFS1CLR = 2; // Clear ADC interrupt flag

  samples[sampleCounter++] = ADC1BUF0 - 512;

  if(sampleCounter >= N_WAVE) {
   sampleCounter = 0;
   samplesReady = true;
   AD1CON1bits.ASAM = 0; // Stop autosampling
  }
 }
}

/* 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
*/

/*
  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 short FIX_MPY(short a, short b)
{
  /* shift right one less bit (i.e. 15-1) */
  int c = ((int)a * (int)b) >> 14;
  /* last bit shifted out = rounding-bit */
  b = c & 0x01;
  /* last shift + rounding bit */
  a = (c >> 1) + b;
  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(short fr[], short fi[], short m, short inverse)
{
  int mr, nn, i, j, l, k, istep, n, scale, shift;
  short 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 =  Sinewave[j+N_WAVE/4];
      wi = -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;
}

Ivan

Wed, 16 Nov 2011 10:22:01 +0000

Hi.

What is the sample rate of this code ? I need an example of high frequency sampling, perhaps 1000 ksps ?

Kind Regards

serveurperso

Wed, 16 Nov 2011 11:07:07 +0000

Hi. What is the sample rate of this code ? I need an example of high frequency sampling, perhaps 1000 ksps ? Kind Regards

Hi,

Device max is 1000ksps, I look 500ksps on the datasheet...?

And I need to learn more on the subject, but all is on the PIC32 datasheet : http://ww1.microchip.com/downloads/en/DeviceDoc/61143E.pdf

Section 22.0 Analog to digital converter

Look formulas "FIGURE 22-3: 500 KSPS SAMPLE RATE CALCULATIONS"

I use the Auto-Sample mode, use formulas with my commented registers values :

// ADC Conversion Clock Select
// TPB = (1 / 80MHz) * 2 = 25ns
// TAD = TPB * (ADCS + 1) * 2 = 25ns * 20 = 500ns
AD1CON3bits.ADCS = 19;

// Auto-Sample Time
// SAMC * TAD = 30 * 500ns = 15uS
AD1CON3bits.SAMC = 30;

I want a #define SAMPLERATE with macros inside the code but I need to learn more on the subject...

Pascal

haz18

Tue, 09 Oct 2012 06:35:51 +0000

Hi there.

Great project! I'm working on something similar too. I noticed that there's no input from your code for the sound analyzer. How did you sample the sound?

majenko

Tue, 09 Oct 2012 10:22:52 +0000

The sampling is interrupt driven.

extern "C" {
 void __ISR(_ADC_VECTOR, ipl6) ADCInterruptHandler() {
  IFS1CLR = 2; // Clear ADC interrupt flag

  samples[sampleCounter++] = ADC1BUF0 - 512;

  if(sampleCounter >= N_WAVE) {
   sampleCounter = 0;
   samplesReady = true;
   AD1CON1bits.ASAM = 0; // Stop autosampling
  }
 }
}

This section of code is triggered every time a single sample is available. That sample is added to the "samples" array, and a new sample is waited for. When the "samples" array is full the "samplesReady" flag is set, and the main program reacts to it and uses the contents of the "samples" array to do the FFT and render the graph.

Of course, it would have been even more efficient to use DMA for the sampling - it could do it entirely in the background then.

haz18

Wed, 10 Oct 2012 06:27:53 +0000

The sampling is interrupt driven.

extern "C" {
void __ISR(_ADC_VECTOR, ipl6) ADCInterruptHandler() {
IFS1CLR = 2; // Clear ADC interrupt flag
samples[sampleCounter++] = ADC1BUF0 - 512;
if(sampleCounter >= N_WAVE) {
sampleCounter = 0;
samplesReady = true;
AD1CON1bits.ASAM = 0; // Stop autosampling
}
}
}

This section of code is triggered every time a single sample is available. That sample is added to the "samples" array, and a new sample is waited for. When the "samples" array is full the "samplesReady" flag is set, and the main program reacts to it and uses the contents of the "samples" array to do the FFT and render the graph. Of course, it would have been even more efficient to use DMA for the sampling - it could do it entirely in the background then.

Yes but where does it input the samples from? I do not see any reference to sample the input from.

and where is the value from this function fix_fft(real, imag, LOG2_N_WAVE, 0); stored in after it has been executed?

majenko

Wed, 10 Oct 2012 08:34:29 +0000

Yes but where does it input the samples from? I do not see any reference to sample the input from.

This bit sets up the ADC and selects the port (RB0) to sample from

// ADC setup
 AD1PCFG = ~1;          // PORTB digital, RB0 analog
 AD1CHSbits.CH0SA = 0;  // Channel 0 positive input is AN0
 //AD1CON2bits.VCFG = 1;  // External VREF+ pin
 AD1CON1bits.SSRC = 7;  // Internal counter ends sampling and starts conversion (auto-convert)
 AD1CON1bits.ASAM = 1;  // Sampling begins when SAMP bit is set

and where is the value from this function fix_fft(real, imag, LOG2_N_WAVE, 0); stored in after it has been executed?

You pass two arrays to the fix_fft function. One of "real" numbers, and one of "imaginary" numbers. Only the "real" array has anything in it, and that is the samples you have taken with the ADC. The fix_fft routine takes those samples, and calculates the FFT real and imaginary components, and places those back in to the real and imaginary arrays.

haz18

Wed, 17 Oct 2012 04:27:14 +0000

Is it possible to modify the codes to make it accept floating-point data for the FFT calculation? I am modifying the existing code to calculate a set of random floating point data (512 values) and output the calculation thru serial comm.

majenko

Wed, 17 Oct 2012 09:23:28 +0000

That fix_fft routine is specifically written to be highly optimised with integer mathematics. While it would be possible to re-write it to use floating point maths it would take a lot of work - many of the operations and techniques are specific to integer maths (shifts, etc) and won't work on floating point.

It would be simpler to either use a routine written specifically for floating point, or convert your floating point values into integer values (fixed point maths) through simple pre-multiplication.

mrkobra79

Wed, 04 Sep 2013 20:21:07 +0000

Hi How to Conect TFT to Chipkit Max32? D0-D15=>Max pin Andreas

mrkobra79

Wed, 04 Sep 2013 22:29:20 +0000

Have a problem with the sketch. Since I do not have the Itdb02_Graph16.h for chipKIT.(now Downloads from Internet) I did the sketch on UTFT to change, got 3 errors.

pic32te.cpp:153:12: error: variable or field '__ISR' declared void pic32te.cpp:153:26: error: 'ip16' was not declared in this scope pic32te.cpp:318:3: error: expected '}' at end of input

serveurperso

Mon, 04 Nov 2013 07:41:03 +0000

Hi, I send you the source and my old lib directory by email.