• State Not Answered
  • Locked Locked
  • Replies 8 replies
  • Subscribers 15 subscribers
  • Views 2588 views
  • Users 0 members are here
Options
Related
How was your experience today?
This discussion has been locked.
You can no longer post new replies to this discussion. If you have a question you can start a new discussion

C51 Cooley-Tukey FFT not giving same output as C89 implementation

Daniel Li

I made a FFT program that takes data from inside a buffer, performs FFT, and spits out the equation form. The C89 (ported from C99) version works well, but when I use the EMF8LB1 board to run it after compiling/linking with C51, the output is different. Why is this happening?

I've compared the codes side by side and as far as I can tell there are no meaningful differences. I can provide .map, .asm, etc. files for the C51 version if that helps.

C89 version:

// Function: _fft
// Uses divide and conquer (Cooley-Tukey) to turn buf (L terms of DTS) into Xk for freqeuncies wK = 2*PI*K/L
void _fft(float bufreal[], float bufimag[], float outreal[], float outimag[], int step)
{
        int i;
        if (step < BUFSIZE) {
                _fft(outreal, outimag, bufreal, bufimag, step * 2);
                _fft(outreal + step, outimag + step, bufreal + step, bufimag + step, step * 2);
                // Do small-size
                for (i = 0; i < BUFSIZE; i += 2 * step) {
                        int j;
                        float treal = outreal[i+step]*cosf(-PI*i/BUFSIZE)-outimag[i+step]*sinf(-PI*i/BUFSIZE);
                        float timag = outimag[i+step]*cosf(-PI*i/BUFSIZE)+outreal[i+step]*sinf(-PI*i/BUFSIZE);
                        printf("\nOut: ");
                        for (j = 0; j < BUFSIZE; j++)
                                if (!outimag[j])
                                        printf("%g ", outreal[j]);
                                else
                                        printf("(%g, %g) ", outreal[j], outimag[j]);
                        printf("\ni: %d; step: %d", i, step);
                        printf("\nout[i+step]: ( %g, %g )", outreal[i+step], outimag[i+step]);
                        printf("\nt: ( %g, %g )", treal, timag);
                        bufreal[i / 2] = outreal[i] + treal;
                        bufimag[i / 2] = outimag[i] + timag;
                        bufreal[(i+BUFSIZE)/2] = outreal[i] - treal;
                        bufimag[(i+BUFSIZE)/2] = outimag[i] - timag;
                }
        }
}

int main()
{
        //PI = atan2(1, 1) * 4; // Define PI
        float bufreal[BUFSIZE] = {1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0};  // Dummy input
        float bufimag[BUFSIZE] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
        float mag[BUFSIZE/2];
        float phs[BUFSIZE/2];
        float outreal[BUFSIZE]; // Create a copy of data
        float outimag[BUFSIZE];
        //showcplx("Data: ", bufreal, bufimag);
        int i;
        float bufsize_dbl = (float) BUFSIZE;
        float sampleLength = (BUFSIZE/SAMPLERATE);

        printf("\rData : ");  // Print string
        // Loop through and print
        for (i = 0; i < BUFSIZE; i++)
                if (!bufimag[i])
                        printf("%g ", bufreal[i]);
                else
                        printf("(%g, %g) ", bufreal[i], bufimag[i]);


        for (i = 0; i < BUFSIZE; i++) {
                outreal[i] = bufreal[i];
                outimag[i] = bufimag[i];
        }
        _fft(bufreal, bufimag, outreal, outimag, 1);

        // magnitude(bufreal, bufimag);
        // Loop over all terms of buf less than folding frequency
        for (i = 0; i < BUFSIZE/2; i++ ) {
                mag[i] = (float) 2/bufsize_dbl*sqrtf( bufreal[i]*bufreal[i] + bufimag[i]*bufimag[i] );
                // Come through too small magnitudes
                if ( mag[i] < THRESH )
                        mag[i] = 0;
        }

        // phase(bufreal, bufimag);
        // Loop over all terms of buf
        for ( i = 0; i < BUFSIZE/2; i++ ) {
                // No magnitude case -> phase = 0
                if ( bufreal[i] < THRESH && bufreal[i] > -THRESH && bufimag[i] < THRESH && bufimag[i] > -THRESH )
                        phs[i] = 0;
                // 90 degree case (real = 0, imag > 0)
                else if ( bufreal[i] < THRESH && bufreal[i] > -THRESH && bufimag[i] > THRESH )
                        phs[i] = 0.5*PI;
                // 180 degree case (real < 0, imag = 0)
                else if ( bufreal[i] < -THRESH && bufimag[i] < THRESH && bufimag[i] > -THRESH )
                        phs[i] = PI;
                // 270 degree case (rea = 0, imag < 0)
                else if ( bufreal[i] < THRESH && bufreal[i] > -THRESH && bufimag[i] < -THRESH )
                        phs[i] = 1.5*PI;
                // All other cases including 0 degree case
                else
                        phs[i] = (float) atan2f(bufimag[i], bufreal[i]);
        }

        // showcplx("\nFFT : ", bufreal, bufimag);
        printf("\nFFT : ");   // Print string
        // Loop through and print
        for (i = 0; i < BUFSIZE; i++)
                if (!bufimag[i])
                        printf("%g ", bufreal[i]);
                else
                        printf("(%g, %g) ", bufreal[i], bufimag[i]);

        // sinusoid();
        printf("\n\n\n"); // Make some space
        // Loop over all terms of mag and phs
        for ( i = 0; i < BUFSIZE/2; i++ ) {
                float omega = 2*PI*i*(1/sampleLength);
                // No magnitude case
                if ( mag[i] == 0.0 )
                        continue; // Do nothing for that term
                // DC term case
                else if ( i == 0 )
                        printf("%f", mag[0]);
                // No phase case
                else if ( phs[i] == 0.0 )
                        printf("+%fcos(%ft)", mag[i], omega);
                // Normal case
                else
                        printf("+%fcos(%ft+%f)", mag[i], omega, phs[i]);
        }

        return 0;
}

C51 version is a carbon copy but with all variables in xdata.

Thanks to everyone in advance!