Dakuan
Anmeldungsdatum: 2. November 2004
Beiträge: 6345
Wohnort: Hamburg
|
Ich habe da ein Programm, das "nebenbei" noch eine große Menge Audiosamples verarbeite soll. Also etwas ähnliches wie ein Spektrogramm. Und das soll irgendwann auch in Echtzeit laufen. Um die Datenmenge, die dabei im Speicher verschoben werden muss, zu halbieren, bin ich auf die Idee gekommen float's zu verwenden. Ist ja eigentlich auch ausreichend. Allerdings wird an einer Stelle im Programm eine vectorielle Addition vorgenommen:
sqrt( re * re + im * im);
Das dabei auftretende Zwischenergebniss passt aber nicht immer in die 23 Bit lange Mantisse eines floats (Eingangswerte: 0 - 32767). Um die Auswirkungen zu testen, habe ich ein Testprogramm erstellt. Das Ergebnis ist, das es kaum nennenswerte Unterschiede gibt.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 | /* Vergleich float vs. double
** Wertebereich: 0 bis 32767
**
** gcc -Wall -o float_test -lm float_test.c
**
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
int main( int argc, char * argv[] ) {
double in1, in2;
double res_d;
float res_f;
float re, im;
if( argc < 2 ) {
printf( "usage: %s value1 value2\n\n", argv[0] );
return 0;
}
in1 = atof( argv[1] );
in2 = atof( argv[2] );
re = (float)in1;
im = (float)in2;
res_d = in1 * in1;
res_f = re * re;
printf( "res_d: %f\nres_f: %f\n\n", res_d, res_f );
res_d = in1 * in1 + in2 * in2;
res_f = re * re + im * im;
printf( "res_d: %f\nres_f: %f\n\n", res_d, res_f );
res_d = sqrt( in1 * in1 + in2 * in2 ); // <<<
res_f = sqrtf( re * re + im * im ); // <<<
printf( "res_d: %f\nres_f: %f\n\n", res_d, res_f );
return 0;
}
|
Habe ich bei meinen Überlegungen einen Denkfehler gemacht oder etwas falsch verstanden? Moderiert von rklm: Wunschgemäß verschoben
|
seahawk1986
Anmeldungsdatum: 27. Oktober 2006
Beiträge: 11180
Wohnort: München
|
Hast du getestet, ob der cast von double zu float auf deiner Zielarchitektur tatsächlich einen Performance-Vorteil gegenüber dem Weiterrechnen mit double-Werten bringt? Wenn man da mal 10000 Paare Zufallszahlen aus dem Wertebereich (mit bis zu 12 Nachkommastellen) über eine abgewandelte Version deines C-Programms jagt, die die Differenz zwischen den mit doubles und floats errechneten Ergebnissen ausgibt, sieht man, dass man bei der letzten Berechnung (Test 2) nur davon ausgehen kann, dass die erste Nachkommastelle bei den beiden Ergebnissen identisch ist.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42 | /* Vergleich float vs. double
** Wertebereich: 0 bis 32767
**
** gcc -Wall -o float_test -lm float_test.c
**
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
int main( int argc, char * argv[] ) {
double in1, in2;
double res_d;
float res_f;
float re, im;
double diff;
if( argc < 2 ) {
printf( "usage: %s value1 value2\n\n", argv[0] );
return 0;
}
in1 = atof( argv[1] );
in2 = atof( argv[2] );
re = (float)in1;
im = (float)in2;
res_d = in1 * in1;
res_f = re * re;
diff = res_d - res_f;
printf("%0.14f;", diff); // Test 0
res_d = in1 * in1 + in2 * in2;
res_f = re * re + im * im;
diff = res_d - res_f;
printf("%0.14f;", diff); // Test 1
res_d = sqrt( in1 * in1 + in2 * in2 ); // <<<
res_f = sqrtf( re * re + im * im ); // <<<
diff = res_d - res_f;
printf("%0.14f\n", diff); // Test 2
return 0;
}
|
Und das Test-Skript in Python3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 | #!/usr/bin/env python3
import random
import statistics
import subprocess
random.seed()
r1 = []
r2 = []
r3 = []
for i in range(10000):
a = random.uniform(0, 32767)
b = random.uniform(0, 32767)
r = subprocess.check_output(['./float_test', str(a), str(b)], universal_newlines=True)
test1, test2, test3 = r.split(";")
r1.append(float(test1))
r2.append(float(test2))
r3.append(float(test3))
for i, t in enumerate((r1, r2, r3)):
print("Test %d| min: %0.14f | max: %0.14f | median: %0.14f" % (i, min(t), max(t), statistics.median(t)))
|
liefert z.B.:
Test 0| min: -91.04566824436188 | max: 91.37555241584778 | median: -0.00029791835550
Test 1| min: -190.70973753929138 | max: 201.20456981658936 | median: -0.00160697591491
Test 2| min: -0.00370822435798 | max: 0.00381295046100 | median: 0.00000391243020
|
eagle87
Anmeldungsdatum: 29. Dezember 2010
Beiträge: 85
|
Hallo Dakuan, ich kann zwar nicht so genau erklären warum, aber ich glaube, dass durch die Quadrierung nur wenige entscheidende Bits im hinteren Teil der Mantisse liegen und abgeschnitten werden. Hab mal deinen Code etwas angepasst um das zu verdeutlichen. 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162 | /* Vergleich float vs. double
** Wertebereich: 0 bis 32767
**
** gcc -Wall -o float_test -lm float_test.c
**
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
void bin_print_float(float val)
{
unsigned char *ptr;
ptr = (unsigned char*) &val;
fputs("\x1B[34m", stdout);
putchar((ptr[3]&0x80)?'1':'0');
fputs(" ", stdout);
fputs("\x1B[32m", stdout);
putchar((ptr[3]&0x40)?'1':'0');
putchar((ptr[3]&0x20)?'1':'0');
putchar((ptr[3]&0x10)?'1':'0');
putchar((ptr[3]&0x08)?'1':'0');
putchar((ptr[3]&0x04)?'1':'0');
putchar((ptr[3]&0x02)?'1':'0');
putchar((ptr[3]&0x01)?'1':'0');
putchar((ptr[2]&0x80)?'1':'0');
fputs("\x1B[31m", stdout);
putchar((ptr[2]&0x40)?'1':'0');
putchar((ptr[2]&0x20)?'1':'0');
putchar((ptr[2]&0x10)?'1':'0');
putchar((ptr[2]&0x08)?'1':'0');
putchar((ptr[2]&0x04)?'1':'0');
putchar((ptr[2]&0x02)?'1':'0');
putchar((ptr[2]&0x01)?'1':'0');
putchar((ptr[1]&0x80)?'1':'0');
putchar((ptr[1]&0x40)?'1':'0');
putchar((ptr[1]&0x20)?'1':'0');
putchar((ptr[1]&0x10)?'1':'0');
putchar((ptr[1]&0x08)?'1':'0');
putchar((ptr[1]&0x04)?'1':'0');
putchar((ptr[1]&0x02)?'1':'0');
putchar((ptr[1]&0x01)?'1':'0');
putchar((ptr[0]&0x80)?'1':'0');
putchar((ptr[0]&0x40)?'1':'0');
putchar((ptr[0]&0x20)?'1':'0');
putchar((ptr[0]&0x10)?'1':'0');
putchar((ptr[0]&0x08)?'1':'0');
putchar((ptr[0]&0x04)?'1':'0');
putchar((ptr[0]&0x02)?'1':'0');
putchar((ptr[0]&0x01)?'1':'0');
fputs("\x1B[0m\n", stdout);
}
void bin_print_double(double val)
{
unsigned char *ptr;
ptr = (unsigned char*) &val;
fputs("\x1B[34m", stdout);
putchar((ptr[7]&0x80)?'1':'0');
fputs("\x1B[32m", stdout);
putchar((ptr[7]&0x40)?'1':'0');
putchar((ptr[7]&0x20)?'1':'0');
putchar((ptr[7]&0x10)?'1':'0');
putchar((ptr[7]&0x08)?'1':'0');
putchar((ptr[7]&0x04)?'1':'0');
putchar((ptr[7]&0x02)?'1':'0');
putchar((ptr[7]&0x01)?'1':'0');
putchar((ptr[6]&0x80)?'1':'0');
putchar((ptr[6]&0x40)?'1':'0');
putchar((ptr[6]&0x20)?'1':'0');
putchar((ptr[6]&0x10)?'1':'0');
fputs("\x1B[31m", stdout);
putchar((ptr[6]&0x08)?'1':'0');
putchar((ptr[6]&0x04)?'1':'0');
putchar((ptr[6]&0x02)?'1':'0');
putchar((ptr[6]&0x01)?'1':'0');
putchar((ptr[5]&0x80)?'1':'0');
putchar((ptr[5]&0x40)?'1':'0');
putchar((ptr[5]&0x20)?'1':'0');
putchar((ptr[5]&0x10)?'1':'0');
putchar((ptr[5]&0x08)?'1':'0');
putchar((ptr[5]&0x04)?'1':'0');
putchar((ptr[5]&0x02)?'1':'0');
putchar((ptr[5]&0x01)?'1':'0');
putchar((ptr[4]&0x80)?'1':'0');
putchar((ptr[4]&0x40)?'1':'0');
putchar((ptr[4]&0x20)?'1':'0');
putchar((ptr[4]&0x10)?'1':'0');
putchar((ptr[4]&0x08)?'1':'0');
putchar((ptr[4]&0x04)?'1':'0');
putchar((ptr[4]&0x02)?'1':'0');
putchar((ptr[4]&0x01)?'1':'0');
putchar((ptr[3]&0x80)?'1':'0');
putchar((ptr[3]&0x40)?'1':'0');
putchar((ptr[3]&0x20)?'1':'0');
putchar((ptr[3]&0x10)?'1':'0');
putchar((ptr[3]&0x08)?'1':'0');
putchar((ptr[3]&0x04)?'1':'0');
putchar((ptr[3]&0x02)?'1':'0');
putchar((ptr[3]&0x01)?'1':'0');
putchar((ptr[2]&0x80)?'1':'0');
putchar((ptr[2]&0x40)?'1':'0');
putchar((ptr[2]&0x20)?'1':'0');
putchar((ptr[2]&0x10)?'1':'0');
putchar((ptr[2]&0x08)?'1':'0');
putchar((ptr[2]&0x04)?'1':'0');
putchar((ptr[2]&0x02)?'1':'0');
putchar((ptr[2]&0x01)?'1':'0');
putchar((ptr[1]&0x80)?'1':'0');
putchar((ptr[1]&0x40)?'1':'0');
putchar((ptr[1]&0x20)?'1':'0');
putchar((ptr[1]&0x10)?'1':'0');
putchar((ptr[1]&0x08)?'1':'0');
putchar((ptr[1]&0x04)?'1':'0');
putchar((ptr[1]&0x02)?'1':'0');
putchar((ptr[1]&0x01)?'1':'0');
putchar((ptr[0]&0x80)?'1':'0');
putchar((ptr[0]&0x40)?'1':'0');
putchar((ptr[0]&0x20)?'1':'0');
putchar((ptr[0]&0x10)?'1':'0');
putchar((ptr[0]&0x08)?'1':'0');
putchar((ptr[0]&0x04)?'1':'0');
putchar((ptr[0]&0x02)?'1':'0');
putchar((ptr[0]&0x01)?'1':'0');
fputs("\x1B[0m\n", stdout);
}
int main( int argc, char * argv[] )
{
double in1, in2;
double res_d;
float res_f;
float re, im;
if( argc < 2 ) {
printf( "usage: %s value1 value2\n\n", argv[0] );
return 0;
}
in1 = atof( argv[1] );
in2 = atof( argv[2] );
re = (float)in1;
im = (float)in2;
res_d = in1 * in1;
res_f = re * re;
bin_print_double(res_d);
bin_print_float( res_f);
printf( "res_d: %f\nres_f: %f\n\n", res_d, res_f );
res_d = in1 * in1 + in2 * in2;
res_f = re * re + im * im;
bin_print_double(res_d);
bin_print_float( res_f);
printf( "res_d: %f\nres_f: %f\n\n", res_d, res_f );
res_d = sqrt( in1 * in1 + in2 * in2 ); // <<<
res_f = sqrtf( re * re + im * im ); // <<<
bin_print_double(res_d);
bin_print_float( res_f);
printf( "res_d: %f\nres_f: %f\n\n", res_d, res_f );
return 0;
}
|
Erst wenn die Wurzel geogen wird gibt es hier größere Unterschiede. Deshalb sind hier auch die Abweichungen etwas größer (von der Anzahl der signifikanten Stellen her betrachtet). Sind die Nachkommastellen für dich relevant? Gruß
Michael
|
Dakuan
(Themenstarter)
Anmeldungsdatum: 2. November 2004
Beiträge: 6345
Wohnort: Hamburg
|
Hast du getestet, ob der cast von double zu float auf deiner Zielarchitektur tatsächlich einen Performance-Vorteil gegenüber dem Weiterrechnen mit double-Werten bringt?
In diese Fall nicht, aber bei meinen Experimenten mit Bildmanipulation. Der erhoffte Performance-Vorteil besteht nicht in der Berechnung der Werte, da man davon ausgehen kann, das die FPU immer mit der vollen Bitbreite rechnet. Es geht viel mehr darum die Ergebnisswerte im Speicher zu verschieben und natürlich auch darum unnötige Konvertierungen zu vermeiden. Also in diesem Fall geht es um die Ergebnisse einer FFT, wo dann das gesamte Ergebnis eines Zeitabschnittes visuell aufbereitet wird. Da aber auch bei der FFT viele Speicherzugriffe erforderlich sind, erhoffe ich mir auch dadurch etwas Zeitvorteil (immerhin halbiert sich auch dabei Anzahl der zu lesenden Bytes).
... sieht man, dass man bei der letzten Berechnung (Test 2) nur davon ausgehen kann, dass die erste Nachkommastelle bei den beiden Ergebnissen identisch ist.
Du willst mir damit sagen, dass ich wieder mal Probleme sehe, wo keine sind? Tatsächlich sind die Nachkommastellen nicht wirklich von Interesse, zumal sich die Eingangsdaten nur im Bereich von 0 bis 32737 bewegen können (16 Bit PCM). Allerdings können bei einigen Filtertypen, die ich einsetze, durch Überschwinger auch höhere Werte auftreten falls meine automatische Pegelbegrenzung wieder mal zu langsam ist. eagle87 schrieb: ... aber ich glaube, dass durch die Quadrierung nur wenige entscheidende Bits im hinteren Teil der Mantisse liegen und abgeschnitten werden.
Sind die Nachkommastellen für dich relevant?
Nicht wirklich. Also pegelmäßig gesehen liegen die unter -90dB. Und für meine Auswertung sind nur Pegel von 20dB unter Maximal Pegel wirklich auswertbar. Um mal zu verdeutlichen, was ich da mache, zeige ich mal im Anhang wie so eine Auswertung aussehen kann (hier noch mit double gerechnet). Jede Pixel Zeile entspricht etwa 1/8 Sek.
|
rleofield
Anmeldungsdatum: 14. September 2008
Beiträge: 779
Wohnort: Görlitz
|
Dakuan schrieb: Ich habe da ein Programm, das "nebenbei" noch eine große Menge Audiosamples verarbeite soll. Also etwas ähnliches wie ein Spektrogramm. Und das soll irgendwann auch in Echtzeit laufen. Um die Datenmenge, die dabei im Speicher verschoben werden muss, zu halbieren, bin ich auf die Idee gekommen float's zu verwenden. Ist ja eigentlich auch ausreichend. Allerdings wird an einer Stelle im Programm eine vectorielle Addition vorgenommen:
sqrt( re * re + im * im);
Das dabei auftretende Zwischenergebniss passt aber nicht immer in die 23 Bit lange Mantisse eines floats (Eingangswerte: 0 - 32767). Um die Auswirkungen zu testen, habe ich ein Testprogramm erstellt. Das Ergebnis ist, das es kaum nennenswerte Unterschiede gibt.
Die Werte kann man auch in 16 bit unsigned aufheben. Nur bei der Rechnung braucht man double. Dass nimmt keinen Platz weg. Wobei man die Multiplikation noch mit 64 bit machen kann, das ist schnell. 1
2
3
4
5
6
7
8
9
10
11
12
13
14 | uint16_t int1 = 32767; // example val 1
uint16_t int2 = 32767; // example val 2
uint64_t res1 = static_cast<uint64_t>(int1);
res1 *= res1;
uint64_t res2 = static_cast<uint64_t>(int2);
res2 *= res2;
uint64_t rr = res1 + res2;
double d = sqrt( static_cast<double>(rr) );
float f = sqrtf( static_cast<float>(rr) );
rr ist 2147352578 = 10 Stellen
f ist 46339.5352 = ca. 7 Stellen genau
d ist 46339.535798279205 = ca. 15 Stellen genau
|
Mal mit dem Debugger anschauen. Der Unterschied liegt im Bereich der möglichen Genauigkeit, mehr ist es nicht.
|
senden9
Anmeldungsdatum: 8. Februar 2010
Beiträge: 965
Wohnort: Österreich
|
Dakuan schrieb: […] Um mal zu verdeutlichen, was ich da mache, zeige ich mal im Anhang wie so eine Auswertung aussehen kann (hier noch mit double gerechnet). Jede Pixel Zeile entspricht etwa 1/8 Sek.
Der Anhang fehlt.
|
Dakuan
(Themenstarter)
Anmeldungsdatum: 2. November 2004
Beiträge: 6345
Wohnort: Hamburg
|
Sorry, aber ich habe es in den letzten 5 Minuten nicht geschafft, den Anhang Hochzuladen, da die Datei angeblich bereits existiert. Letzter Versuch ...
- Bilder
|
Vain
Anmeldungsdatum: 12. April 2008
Beiträge: 2503
|
(Möchte nur gerade semi-OT in den Raum werfen, dass fftw3 eine okay benutzbare Bibliothek für FFTs ist. Möglicherweise ist das eine Option für dich, sodass du das gar nicht selbst schreiben musst.)
|
Dakuan
(Themenstarter)
Anmeldungsdatum: 2. November 2004
Beiträge: 6345
Wohnort: Hamburg
|
eagle87 schrieb: ich kann zwar nicht so genau erklären warum, aber ich glaube, dass durch die Quadrierung nur wenige entscheidende Bits im hinteren Teil der Mantisse liegen und abgeschnitten werden.
Ich habe nochmal darüber nachgedacht und denke dass bei dem angenommenen maximalen Eingangswert von 32767 8 Bit verloren gehen (eigentlich 9, aber die führende 1 wird ja nicht mit gespeichert). Dann müsste der maximale Fehler bei
sqrt(255) = 16
liegen, tut er aber offenbar nicht. rleofield schrieb: Die Werte kann man auch in 16 bit unsigned aufheben. Nur bei der Rechnung braucht man double.
Die Werte sind ja ursprünglich 16 Bit Werte. Aber das Programm macht ja noch mehr damit. Die Werte werden z.B. auch gefiltert und mit Sinus und Cosinus multipliziert u.s.w. Da müsste dann doch wieder umgewandelt werden. Die Werte werden von zentraler Stelle an mehrere Module verteilt und daher möchte ich möglichst durchgängig mit einem einheitlichen Format arbeiten. Da bleibt dann wohl nur entweder darauf hoffen, das die Abweichung sich nicht auswirkt oder die Werte nur für diese eine Zeile in double und zurück konvertieren.
|
Dakuan
(Themenstarter)
Anmeldungsdatum: 2. November 2004
Beiträge: 6345
Wohnort: Hamburg
|
@Vain
Auf FFTW habe ich bewusst verzichtet. Und der Standard Code für eine FFT ist auch nur zwischen 50 und 60 Zeilen lang (den habe ich bei Dr.Dobb's abgeschrieben). Das Problem tritt an einer anderen Stelle auf. Es handelt sich da um 2 Datenströme, die um 90 Grad gegeneinander phasenverschoben sind, und wo am Ende der Betrag gebildet werden soll.
|
Dakuan
(Themenstarter)
Anmeldungsdatum: 2. November 2004
Beiträge: 6345
Wohnort: Hamburg
|
@Vain Aufgrund Deines Einwands habe ich mir noch einmal angesehen was da nach der FFT passiert. Man sollte nicht einfach Code Teile aus anderen Projeken recyclen, ohne zu prüfen ob das Umfeld noch stimmt. Tatsächlich tritt hier das selbe Problem auf. Das komplexe Ergebnis der FFT wird auf die gleiche Weise in Betragswerte umgewandelt. Dass es da auch ein Problem geben könnte ist mir allerdings bisher nicht aufgefallen, da die Ergebnisse anschließend normalisiert und für die Anzeige noch einmal skaliert werden (beides nicht Bestandteil der FFT). Ich muss jetzt erstmal austesten, ob in wieweit das einen erkennbaren Einfluss auf das Ergebnis hat und wie groß oder klein der Unterschied bei der Verwendung von floats ist. Meine Großen Rechner merken den Unterschied wahrscheinlich nicht, aber mein Raspi könnte da schon ins schwitzen kommen.
|