4
\$\begingroup\$

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;
ADIE=0;
// 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++)
    {
        // Perform the ADC conversion
        // Select the desired ADC and start the conversion
        ADCON0 = 0b00000011;    // Start the ADC conversion on AN0

        // 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)
        //realNumbers[i] = ((short)(ADRESH << 8) + (short)ADRESL) - 512;
        realNumbers[i] = ((short)(ADRESH << 8) + (short)ADRESL)-2048;

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

\$\endgroup\$
  • \$\begingroup\$ 201 was generated like this: table[index] = (sin(index * 90/(N_WAVE/4)) * 32768). table[1] = 201 \$\endgroup\$ – kkrambo Oct 23 '15 at 14:10
  • \$\begingroup\$ 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. \$\endgroup\$ – kkrambo Oct 23 '15 at 14:12
  • \$\begingroup\$ Thanks. I'll see if recreating the lookup table solves my problem. \$\endgroup\$ – cheunste Oct 24 '15 at 3:38
1
\$\begingroup\$

You're confusing number of samples with resolution.

The 16-bit transform means each sample is 16-bits. In other words, the samples can be -32768 to 32767. The number of samples (N_WAVE) taken is completely irrelevant to the 'bitness' of the FFT. The 16-bits is referring to sample size, not the number of samples being taken. So the 1024 sample count and the resolution of those said samples have no relationship to each other, and 'where it came from' is, well, from the number of samples you wanted to take of a wave. This is determined by the nyquist frequency. You can just as readily use an N_WAVE of 256 if you wanted.

As for why your code isn't working, you're FFT is probably overflowing like crazy. The PIC doesn't magically change multiplication of 2 16-bit values into a 32 bit value like it does for turning 2 8-bit values into a 16-bit value. You need more than 16-bits of storage if you're multiplying 16-bit numbers. Even if your dynamic range is just 12 bits, you're scaling it, and if you didn't, well, anything except 0 multiplied by some of the values in your lookup table (which span the entire signed short value range) will overflow.

In the fix_fft function, use ints instead of shorts for the two arrays and m. This will force 32 bit multiplication and will resolve any overflow issues.

Honestly, this is probably not a question suited to electronics, so you might get a better answer in a programming related stackexchange. But I think this will help. Good luck!

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.