Mail Archives: djgpp/2004/03/02/20:27:50
Dear Hans and Listmembers -
>
>We haven't heard any feedback from you in ages... did you try any
>of the suggestions? Did they work?
>[...]
>
>--
>Hans-Bernhard Broeker (broeker AT physik DOT rwth-aachen DOT de)
>Even if all the snow were burnt, ashes would remain.
Please pardon the delay. I appreciate EVER SO MUCH your suggestions. To
the best of my understanding, I have tried to implement them. The answer
to the problem is still elusive...
I have compiled the program using double and long double. The results follow:
[note HBR is my DOS computer; SARAH is my wife's Win2000 computer]
When the programs were run, on the windows OS, nothing else was running on
windows. I have yet to run the programs on DOS 7.0 [the pure DOS shell of
Win 98]
[min is time, delta is the calculated value of interest, etot is the
integral under the FFT curve]
HBR/float
min delta etot
0.45 141.79744 490616.065924
3.15 -3.106526 1270.050473
SARAH/float
min delta etot
0.45 5.886202 56.013115
3.15 -3.106526 1270.050473
HBR/double
min delta etot
0.45 102.283497 895621.961731
3.15 -3.106527 1270.050546
SARAH/double
min delta etot
0.45 0.399665 187.450401
3.15 -3.106527 1270.050546
HBR/long double
min delta etot
0.45 135.183805 650175.488059
3.15 -3.106527 1270.050546
SARAH/long double
min delta etot
0.45 135.183805 650175.488059
3.15 -3.106527 1270.050546
The problem with the long double calculations, I think, is that the sqrt
function assumes double for input.
The FFT routine which I have not modified at all is:
/*
* Here is the real work-horse.
* It's a generic FFT, so nothing is lost or approximated.
* The samples in wave[] should be in order, and they
* will be decimated when fft() returns.
*/
void fft (complx wave[], int bits)
{
register int loop, loop1, loop2;
unsigned i1; /* going to right shift this */
int i2, i3, i4, y;
FLOAT a1, a2, b1, b2, z1, z2;
i1 = FFT_LEN / 2;
i2 = 1;
/* perform the butterfly's */
for (loop = 0; loop < bits; loop++)
{
i3 = 0;
i4 = i1;
for (loop1 = 0; loop1 < i2; loop1++)
{
y = PERMUTE(i3 / (int)i1, bits);
z1 = CosineTable[y];
z2 = -SineTable[y];
for (loop2 = i3; loop2 < i4; loop2++)
{
a1 = REAL(loop2);
a2 = IMAG(loop2);
b1 = z1 * REAL(loop2+i1) - z2 * IMAG(loop2+i1);
b2 = z2 * REAL(loop2+i1) + z1 * IMAG(loop2+i1);
REAL(loop2) = a1 + b1;
IMAG(loop2) = a2 + b2;
REAL(loop2+i1) = a1 - b1;
IMAG(loop2+i1) = a2 - b2;
}
i3 += (i1 << 1);
i4 += (i1 << 1);
}
i1 >>= 1;
i2 <<= 1;
}
}
Recognizing the the problem may be from overflow of the float, certain
changes were made in the program:
The original amplitude calculating functions:
/*
* 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)); */
}
Changed to:
/*
* 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*(double)a) + ((double)b*(double)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 (hypot (REAL(n), IMAG(n)));
/* return (REAL(n), IMAG(n)); */
}
>>>NOTE the hypot2 was ONLY used for float calls <<<<<
>>> If the problem is in the c*c = a*a + b*b, then I could use a routine
for that in Long Double
<<<<
The original integration routine:
/* 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;
}
Changed to:
[to sort the numbers prior to the sum]
/* 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, temp3[1450];
int i, j, jj;
sum = 0.0;
for (i = 0; i < 1450; i++)
temp3[i] = 0.0;
sum += temp2[start_pt];
/* for(i = start_pt+1; i < end_pt; i++)
sum += 2*temp2[i];*/
jj = 0;
for(i = start_pt+1, j = 0; i < end_pt; i++, j++)
{
temp3[j] = temp2[i];
jj = j;
}
qsort(temp3, jj, sizeof(FLOAT), mycomp);
for(i = 0; i < jj; i++)
sum += 2*temp3[i];
sum += temp2[end_pt];
sum = (end_val-start_val)/(2*((FLOAT)end_pt-(FLOAT)start_pt))*sum;
return sum;
}
=============
I read the paper "What every computer Scientist needs to know about
floating point..". My background is NOT in numerical analysis, so, I do
not understand it totally. I hope I can find a mathematician that can help.
MUCH THANKIS IN ADVANCE FOR ALL YOUR EXCELLENT HELP.
Ethan Rosenberg
- Raw text -