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

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!

0 Daniel Li over 6 years ago

C51 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 xdata i;
        //int xdata j;
        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 xdata j;
                        float xdata treal = outreal[i+step]*cosf(-PI*i/BUFSIZE)-outimag[i+step]*sinf(-PI*i/BUFSIZE);
                        float xdata 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[i], outimag[i]);
                        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 xdata bufreal[BUFSIZE] = {1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0};    // Dummy input
        float xdata bufimag[BUFSIZE] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
        float xdata mag[BUFSIZE/2];
        float xdata phs[BUFSIZE/2];
        float xdata outreal[BUFSIZE];   // Create a copy of data
        float xdata outimag[BUFSIZE];
        //showcplx("Data: ", bufreal, bufimag);
        int xdata i;
        float xdata bufsize_dbl = (float) BUFSIZE;
        float xdata 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 xdata 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;
}

<\pre>

0 Daniel Li over 6 years ago in reply to Daniel Li

C89 version output (verified to be correct)

Data : 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
Out: 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
i: 0; step: 8
out[i+step]: ( 1, 0 )
t: ( 1, 0 )
Out: 1 0 -1 0 1 0 -1 0 1 0 -1 0 (5.89021e-039, 1) 2 (16, -1) 2.24208e-044
i: 0; step: 8
out[i+step]: ( 1, 0 )
t: ( 1, 0 )
Out: 2 0 -1 0 2 0 -1 0 0 0 -1 0 0 0 -1 0
i: 0; step: 4
out[i+step]: ( 2, 0 )
t: ( 2, 0 )
Out: 2 0 -1 0 2 0 -1 0 0 0 -1 0 0 0 -1 0
i: 8; step: 4
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 (5.89021e-039, 4) 2
i: 0; step: 8
out[i+step]: ( -1, 0 )
t: ( -1, 0 )
Out: -1 0 0 0 -1 0 0 0 -1 0 (5.89021e-039, 4) 2 (16, -1) 2.24208e-044 8.99954e-039 5.88424e-039
i: 0; step: 8
out[i+step]: ( -1, 0 )
t: ( -1, 0 )
Out: -2 0 2 0 -2 0 0 0 0 0 0 0 0 0 (1.01064e-038, 2) 2.32438e-017
i: 0; step: 4
out[i+step]: ( -2, 0 )
t: ( -2, 0 )
Out: -2 0 2 0 -2 0 0 0 0 0 0 0 0 0 (1.01064e-038, 2) 2.32438e-017
i: 8; step: 4
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 4 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0
i: 0; step: 2
out[i+step]: ( -4, 0 )
t: ( -4, 0 )
Out: 4 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0
i: 4; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 4 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0
i: 8; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 4 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0
i: 12; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4)
i: 0; step: 8
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4) 2 (16, -4) 2.24208e-044 8.99954e-039
i: 0; step: 8
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 1.01064e-038
i: 0; step: 4
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 1.01064e-038
i: 8; step: 4
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 0 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4) 2 (16, -4)
i: 0; step: 8
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4) 2 (16, -4) 2.24208e-044 8.99954e-039 5.88424e-039 8.99963e-039
i: 0; step: 8
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 8 0 0 0 0 0 0 0 1.01064e-038 2.32438e-017 8.9994e-039
i: 0; step: 4
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 8 0 0 0 0 0 0 0 1.01064e-038 2.32438e-017 8.9994e-039
i: 8; step: 4
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4)
i: 0; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4)
i: 4; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4)
i: 8; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 -4 0 0 0 0 0 0 0 0 0 0 0 0 0 (5.89021e-039, 4)
i: 12; step: 2
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 0; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 2; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 4; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 6; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, 0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 8; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 10; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 12; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
Out: 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0
i: 14; step: 1
out[i+step]: ( 0, 0 )
t: ( 0, -0 )
FFT : 0 0 0 0 8 0 0 0 0 0 0 0 8 0 0 0


+1.000000cos(12.566371t)
Press any key to continue . . .

0 Daniel Li over 6 years ago in reply to Daniel Li

C51 output (not the same values for some reason):

Data : 1.000000 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00
Out: 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
i: 0; step: 32
out[i+step]: ( 1.000000, 0.000000e+00 )
t: ( 1.000000, 0.000000e+00 )
Out: 1.000000 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00
i: 0; step: 64
out[i+step]: ( 0.000000e+00, 0.000000e+00 )
t: ( 0.000000e+00, 0.000000e+00 )
Out: 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
i: 0; step: 128
out[i+step]: ( 0.000000e+00, 0.000000e+00 )
t: ( 0.000000e+00, 0.000000e+00 )
Out: 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
i: 0; step: 256
out[i+step]: ( 0.000000e+00, 5.911820e-39 )
t: ( 0.000000e+00, 0.000000e+00 )
FFT : (1.000000, 1.000000) 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00 (1.000000, -1.000000) 0.000000e+00 -1.000000 0.000000e+00 1.000000 0.000000e+00 -1.000000 0.000000e+00


0.176776+0.125000cos(6.283185t+3.141592)+0.125000cos(12.566371t)+0.125000cos(18.849554t+3.141592)

0 HansBernhard Broeker over 6 years ago

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

You'll want to double-check that fact. Get a file comparison tool, and use it.
Cancel
Vote up 0 Vote down

Cancel
0 Ima Spammer over 6 years ago in reply to HansBernhard Broeker

Remember that Keil C51 has many "tricks" to make standard 'C' work on the weird and restrictive 8051 architecture.

So make sure you are in "pure ANSI" mode - and not using any of the special options such as disabling integer promotion ...

Apart from that, just run up 2 instances of uVision and step through the 2 code versions side-by-side in the simulator ...
Cancel
Vote up 0 Vote down

Cancel
0 ²erik malund over 6 years ago in reply to Ima Spammer

C51 version is a carbon copy but with all variables in xdata.
by large model or by specifying at variables?

if by large model, you have DRAMATICALLY changed the timing, could it be an ISR ?

BTW why ?
Cancel
Vote up 0 Vote down

Cancel
0 Daniel Li over 6 years ago in reply to ²erik malund

I did it at each variable i.e. char xdata hello_world.

Thanks everyone for the responses, I ended up replacing the recursive algorithm with a brute force iterative one, since there was not enough ram on my EFM8 to hold the stack required by the recursive algorithm. Now everything works.
Cancel
Vote up 0 Vote down

Cancel
0 Ima Spammer over 6 years ago in reply to Daniel Li

You do realise that you can't just do recursion in C51 as in Standard 'C', don't you?

That's one of the "tricks" I referred to earlier.

http://www.keil.com/support/man/docs/c51/c51_le_reentrantfuncs.htm
Cancel
Vote up 0 Vote down

Cancel