¿La tabla de búsqueda en este código FFT está causando problemas con mi ADC de 12 bits?

Intento implementar un proyecto de analizador de espectro que encontré en aquí , pero en lugar de usar el microcontrolador PIC18F4550, estoy usando un PIC18F46K80.

Al comparar estos dos chips, una diferencia importante que encontré es que el 4550 usa un ADC de 10 bits, mientras que el 46K80 I'm usa el ADC de 12 bits. Sin embargo, me estoy topando con un problema en el que los resultados de la función FFT me dan una matriz con 215, o algunos otros valores aleatorios para todos los grupos de frecuencias y no solo el grupo de 1kHz (vea el código para lo que quiero decir).

Por lo que sé, la forma en que funciona este código es mediante

1) el microcontrolador lee 64 muestras cada 50 us (para obtener una frecuencia de muestreo de 20 kHz) y almacena los resultados en una matriz (números reales).

2) realiza una FFT de 16 bits en la matriz realNumbers y almacena el resultado de la FFT en una matriz pero organizada por grupos de frecuencia.

3) Las matrices reales y complejas se utilizan para calcular la magnitud y se almacenan en una matriz.

main.c Tenga en cuenta que los NOP () solo están allí para propósitos de depuración     // FOSC = HS2     vacío principal     {     // Configuración del puerto PIC -------------------------------------------- ------------     // Configurar el ADC a bordo     // Vss y Vdd como referencias de voltaje

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     #incluir     #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 (donde está contenida la tabla de búsqueda de seno)

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

Al profundizar en el código FFT, solo me estoy descarrilando, pero estoy un poco convencido de que la tabla de búsqueda de onda sinusoidal me está causando problemas, por lo que me viene a la mente lo siguiente:

1) Para una FFT de 16 bits, ¿por qué la tabla de búsqueda Sinewave utiliza una N_WAVE de 10 bits (1024)? ¿De dónde vino esto de repente?

2) ¿Puede la tabla de búsqueda Sinewave darme problemas cuando estoy usando ADC de 12 bits como entrada? Ya he intentado revertir la generación de números en la tabla de búsqueda de onda sinusoidal dejando que sin (pi / 2) = 1 = 32768 y sin (pi / 4) = 0.707 = 23166 pero todavía no puedo terminar con cómo era un 201 generado

3) ¿Crees que es la tabla de búsqueda un posible problema o podría ser el tamaño de muestreo de 64?

Por favor, hágame saber si algo necesita una aclaración.