# Is the lookup table in this FFT code causing problems with my 12-bit ADC?

I'm trying to implement a spectrum analyzer project I found from here, but instead of using the PIC18F4550 microcontroller, I'm using a PIC18F46K80.

When comparing these two chips, a major relevant difference I found is that the 4550 uses a 10-bit ADC while the 46K80 I'm uses 12-bit ADC. However, I'm running into an issue where the results from the FFT function is giving me an array with 215, or some other random values for all frequency buckets and not just the 1kHz bucket (see code for what I mean by this).

As far as I know, the way this code works is by

1) microcontroller reads 64 samples every 50us (to get 20kHz sampling rate) and stores results into an array (realNumbers).

2) performs a 16-bit FFT on the realNumbers array and stores the FFT result in an array but organized by frequency buckets.

3) Both the real and complex arrays are then used to calculate the magnitude and stored in an array.

main.c Note the NOP() are only there for debugging purposes //FOSC=HS2 void main(void) { // PIC port set up -------------------------------------------------------- // Configure on-board ADC // Vss and Vdd as voltage references

ADCON1 = 0b00000000;
// Configure the ADC acquisition time according to the datasheet
ADCON2 = 0b10110110; // Note: output is right justified

// Configure ports as inputs (1) or outputs(0)
TRISA = 0b00000001;
TRISB = 0b00000000;
TRISC = 0b00000011;
TRISD = 0b00000000;
TRISE = 0b00000000;
while(1)
{
// Perform the FFT

// Get 64 samples at 50uS intervals
// 50uS means our sampling rate is 20KHz which gives us
// Nyquist limit of 10Khz
short i = 0;
unsigned short result;
for (i = 0; i < 64; i++)
{
// Select the desired ADC and start the conversion

// Wait for the ADC conversion to complete
TESTPIN_W4 = 1; // Don't remove this... it will affect the sample timing
while(GODONE);
TESTPIN_W4 = 0; // Don't remove this... it will affect the sample timing

// Get the 10-bit ADC result and adjust the virtual ground of 2.5V
// back to 0Vs to make the input wave centered correctly around 0
// (i.e. -512 to +512)

// Set the imaginary number to zero
imaginaryNumbers[i] = 0;

// This delay is calibrated using an oscilloscope according to the
// output on RA1 to ensure that the sampling periods are correct
// given the overhead of the rest of the code and the ADC sampling
// time.
//
// If you change anything in this loop or use the professional
// (optimised) version of Hi-Tech PICC18, you will need to re-
// calibrate this to achieve an accurate sampling rate.
__delay_us(7);
}
// Perform the (forward) FFT calculation

// Note: the FFT result length is half of the input data length
// so if we put 64 samples in we get 32 buckets out.  The first bucket
// cannot be used so in reality our result is 31 buckets.
//
// The maximum frequency we can measure is half of the sampling rate
// so if we sample at 20Khz our maximum is 10Khz (this is called the
// Nyquist frequency).  So if we have 32 buckets divided over 10Khz,
// each bucket represents 312.5Khz of range, so our lowest bucket is
// (bucket 1) 312.5Hz - 625Hz and so on up to our 32nd bucket which
// is 9687.5Hz - 10,000Hz

//  1 : 312.5 - 625
//  2 : 625 - 937.5
//  3 : 937.5 - 1250
//  4 : 1250 - 1562.5
//  5 : 1562.5 - 1875
// .....
// 30 : 9375 - 9687.5
// 31 : 9687.5 - 10000

// Note: the '6' is the size of the input data (2 to the power of 6 = 64)
NOP();

TESTPIN_W5 = 1;
fix_fft(realNumbers, imaginaryNumbers, 10);
NOP();
// Take the absolute value of the FFT results

// Note: The FFT routine returns 'complex' numbers which consist of
// both a real and imaginary part.  To find the 'absolute' value
// of the returned data we have to calculate the complex number's
// distance from zero which requires a little pythagoras and therefore
// a square-root calculation too.  Since the PIC has multiplication
// hardware, only the square-root needs to be optimised.

long place, root;
unsigned int z;
for (char k=1; k < 32; k++)
{
z = (realNumbers[k] * realNumbers[k]);
z = z + (imaginaryNumbers[k] * imaginaryNumbers[k]);

place = 0x40000000;
root = 0;
while (place > z) place = place >> 2;

while (place)
{
if (z >= root + place)
{
z -= root + place;
root += place * 2;
}
root = root >> 1;
place = place >> 2;
}
realNumbers[k] = root;
}
TESTPIN_W5 = 0;

// Now we have 32 buckets of audio frequency data represented as
// linear intensity in realNumbers[]
//
// Since the maximum input value (in theory) to the SQRT function is
// 32767, the peak output at this stage is SQRT(32767) = 181.

// Draw a bar graph of the FFT output data
TESTPIN_W6 = 1;
NOP();
//drawFftGraph(realNumbers);
TESTPIN_W6 = 0;
}
}


fft.c #ifndef FFT_C #define FFT_C #include #include "fft.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.
//
/*
fix_fft() - perform forward fast Fourier transform.
fr[n],fi[n] are real and imaginary arrays, both INPUT AND
RESULT (in-place FFT), with 0 <= n < 2**m
*/
void fix_fft(short fr[], short fi[], short m)
{
long int mr = 0, nn, i, j, l, k, istep, n, shift;
short qr, qi, tr, ti, wr, wi;

n = 1 << m;
nn = n - 1;

/* max FFT size = N_WAVE */
//if (n > N_WAVE) return -1;

/* 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)
{
/*
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.

It may not be obvious, but the shift will be
performed on each data point exactly once,
during this pass.
*/

// Variables for multiplication code
long int c;
short b;

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];

wr >>= 1;
wi >>= 1;

for (i=m; i<n; i+=istep)
{
j = i + l;

// Here I unrolled the multiplications to prevent overhead
// for procedural calls (we don't need to be clever about
// the actual multiplications since the pic has an onboard
// 8x8 multiplier in the ALU):

// tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
c = ((long int)wr * (long int)fr[j]);
c = c >> 14;
b = c & 0x01;
tr = (c >> 1) + b;

c = ((long int)wi * (long int)fi[j]);
c = c >> 14;
b = c & 0x01;
tr = tr - ((c >> 1) + b);

// ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
c = ((long int)wr * (long int)fi[j]);
c = c >> 14;
b = c & 0x01;
ti = (c >> 1) + b;

c = ((long int)wi * (long int)fr[j]);
c = c >> 14;
b = c & 0x01;
ti = ti + ((c >> 1) + b);

qr = fr[i];
qi = fi[i];
qr >>= 1;
qi >>= 1;

fr[j] = qr - tr;
fi[j] = qi - ti;
fr[i] = qr + tr;
fi[i] = qi + ti;
}
}

--k;
l = istep;
}
}

#endif


fft.h (Where the sine lookup table is contained)

#ifndef _FFT_H
#define _FFT_H

// Definitions
#define N_WAVE      1024    // full length of Sinewave[]
#define LOG2_N_WAVE 10      // 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 short Sinewave[N_WAVE-N_WAVE/4] = {
0,    201,    402,    603,    804,   1005,   1206,   1406,
1607,   1808,   2009,   2209,   2410,   2610,   2811,   3011,
3211,   3411,   3611,   3811,   4011,   4210,   4409,   4608,
4807,   5006,   5205,   5403,   5601,   5799,   5997,   6195,
//.....
-32412, -32441, -32468, -32495, -32520, -32544, -32567, -32588,
-32609, -32628, -32646, -32662, -32678, -32692, -32705, -32717,
-32727, -32736, -32744, -32751, -32757, -32761, -32764, -32766,
};

// Function prototypes
void fix_fft(short fr[], short fi[], short m);


Looking deep into the FFT code, I'm only getting more derailed, but I'm kinda convinced that the sinewave lookup table is causing problems for me, which is why the following comes to mind:

1) For a 16-bit FFT, why is the Sinewave lookup table using a N_WAVE of 10 bit (1024)? Where did this suddenly come from?

2) Can the Sinewave lookup table give me issues when I'm using 12-bit ADC as an input? I have already tried to reverse generate the numbers in the sinewave lookup table by letting sin(pi/2)=1=32768 and sin(pi/4)=0.707 = 23166 but I still can't end up with how a 201 was generated

3) Do you think it is the lookup table that's a possible issue or could it be the 64 sampling size?

Please let me know if something needs clarification.

• 201 was generated like this: table[index] = (sin(index * 90/(N_WAVE/4)) * 32768). table[1] = 201 – kkrambo Oct 23 '15 at 14:10
• That's 90 degrees, which is 1/4 of a full cycle. If you're calculating the sine in radians then you'd have to replace 90 degrees with pi/2 radians. – kkrambo Oct 23 '15 at 14:12
• Thanks. I'll see if recreating the lookup table solves my problem. – cheunste Oct 24 '15 at 3:38