Mail Archives: djgpp/2004/02/16/23:13:55
Dear Listmembers -
MUCH THANKS IN ADVANCE FOR ANY AVICE WHICH WILL BE GIVEN.
>Does that mean that the modified program produces _exactly_ the same
>results as before, on each platform? Or do you see any changes in the
>results with the _CRT0_FLAG_FILL_DEADBEEF bit set?
No change wantsoever!!!
The following are code FRAGMENTS from the program. My feelings are that
the error is in the integration routine. etot is defined as FLOAT.
I have extracted these fragments for a program of approx. 3000 lines of
code. I hope I am not burdening you with too much.
************
Ethan Rosenberg
*********************
/* fft.h: Defs for fft routines
* [Part of simple-fft-1.5.tar.Z (ver 1.5)]*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <conio.h>
#include <stdlib.h>
#include <pc.h>
/*
* If you want double-precision floats used everywhere,
* define this to be "double"
*/
#define FLOAT float
/*
* Since every math.h seems to have their own idea of what the
* elements in struct complex are called, we define our own.
*/
struct _complx
{
FLOAT re;
FLOAT im;
};
typedef struct _complx complx;
/*
* This PERMUTE could be defined as a table lookup (Sampson did this),
* but I found that to be significantly slower than just computing it
* inline. Of course, this is w/ GCC, and not MS-DOS Turbo C.
* I'll leave it as a macro; you can tweak with it.
*/
#define PERMUTE(x, y) reverse((x), (y))
#define REAL(x) wave[(x)].re
#define IMAG(x) wave[(x)].im
#define ALPHA 0.54 /* ala hanning */
#define FFT_LEN 16384
/*
* Bit reverser for unsigned ints
* Reverses 'bits' bits.
*/
inline const unsigned int reverse (unsigned int val, int bits)
{
unsigned int retn = 0;
while (bits--)
{
retn <<= 1;
retn |= (val & 1);
val >>= 1;
}
return (retn);
}
/*
* Calculate the hypotenuse of a right triangle. The function in math.h
* assumes the the input and return values are double. This uses FLOAT
*/
FLOAT hypot2(FLOAT a, FLOAT b)
{
FLOAT c,asq,bsq;
double aa,bb,cc;
asq = a*a;
bsq = b*b;
aa = (double)asq;
bb = (double)bsq;
cc = sqrt(aa+bb);
c = (FLOAT)sqrt((double)(a*a) + (double)(b*b));
return c;
}
/*
* Calculate amplitude of component n in the decimated wave[] array.
*/
FLOAT amp (int n, complx wave[], int bits)
{
n = PERMUTE (n, bits);
return (hypot2 (REAL(n), IMAG(n)));
/* return (REAL(n), IMAG(n)); */
}
/* Integrate the PSD function, using the trapaziodal rule */
FLOAT trap_integ(FLOAT temp2[],int start_pt,int end_pt,FLOAT
start_val,FLOAT end_val)
{
FLOAT sum;
int i;
sum = 0.0;
sum += temp2[start_pt];
for(i = start_pt+1; i < end_pt; i++)
sum += 2*temp2[i];
sum += temp2[end_pt];
sum = (end_val-start_val)/(2*((FLOAT)end_pt-(FLOAT)start_pt))*sum;
return sum;
}
...snip ....
for (l = 0; l < FFT_LEN; l++)
fft_out[l] = amp(l, wave, 14);
for(l = 0; l < FFT_LEN/2; l++)
temp2[l] += fft_out[l]*fft_out[l]*2;
e_tot = trap_integ(temp2,103,1537,0.2,3.0);
- Raw text -