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!
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>
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 . . .
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)