Forum: Digitale Signalverarbeitung / DSP / Machine Learning FFT: Numerische Genauigkeit?


von Theodor (Gast)


Lesenswert?

Hi, wenn ich eine FFT mit 2^x Stützstelen nach dem Radix-2-Alogrithmus 
auf Basis von 16-Bit-Zahlen berechnen lasse, wie kann ich dann die 
numerische
Genauigkeit berechnen? Anders: wieviele Bits im Ergebnis kann ich 
wegschmeißen?

von Joe J. (j_955)


Lesenswert?

Hi,

genau darüber zerbreche ich mir auch den Kopf. Die Frage ist, wie du die 
FFT implementierst(Sprache,Plattform). Ich bin dabei, das ganze in C# zu 
machen, und dies ist wesentlich bequemer, weil dort Werkzeuge vorhanden 
sind, welche z.B. die Rundungsfehler dokumentieren.(zumindest im VS)

Grüße

von Purzel H. (hacky)


Lesenswert?

>Die Frage ist, wie du die FFT implementierst (Sprache,Plattform).

Eher nicht. Zuerst muss man die Theorie begriffen haben. Dann kann man 
voraussagen welche Fehler wie propagieren.
Alternativ kann man duch Tests herausfinden wie ein 2bit Rauschen das 
Spektrum veraendert.

von Detlef _. (detlef_a)


Lesenswert?

Für die FFT muß Du die (16 Bit-) Eingangswerte ja mit cos/sin 
Koeffizienten multiplizieren. Die sind in Integer-Darstellung (aber auch 
in float-Darstellung) ja notwendigerweise quantisiert, üblich sind 
Coeffizientenbreiten, die die Hardware verarbeiten kann, also z.B. 16Bit 
oder 18Bit. Das Ergebnis wird dann wieder runterskaliert, also die 
unteren Bits werden weggeschmissen.

Eine 2^x - sample FFT hat x 'stages', die jeweils die Multiplikation mit 
sin/cos machen und eine Addition. Falls die Multiplikation wieder auf 
die ursprüngliche Bitbreite skaliert kann trotzdem pro Stage durch die 
Addition ein Bit 'hinzukommen'. Dann besteht die Gefahr des Überlaufs. 
Also entweder mit höheren Bitbreiten rechnen oder nach jeder stage 
abhängig von den Ergebnissen skalieren. Dieser Bitgewinn reflektiert, 
dass sich zB für einen Sinus die Energie in einer Spektrallinie 
konzentriert.

Also: die Koeffizienten geeignet skalieren, pro stage ein Bit 'Gewinn' 
einrechnen, bißchen Reserve lassen, ggf. dynamisch skalieren. Bei 16Bit 
Eingangsdaten beträgt Dein SNR sowieso maximal ca. 96dB, da machts dann 
keinen Sinn mit 32 oder sogar 40 Bit Zahlen zu rechnen.

Cheers
Detlef

von Andreas S. (andreas) (Admin) Benutzerseite


Lesenswert?

Jo-hann Ha-nnes schrieb:
> Ich bin dabei, das ganze in C# zu machen

Da kannst du einfach float oder double verwenden und brauchst dir keine 
Gedanken mehr über die Genauigkeit machen.

Ein Trick den man bei Festkommazahlen anwenden kann ist, die Daten vor 
der FFT durch Schieben so zu skalieren dass der Wertebereich voll 
ausgenutzt wird. Auf diese Weise kann man auch bei Signalen mit großem 
Dynamikbereich oft mit relativ kurzen Wortlängen auskommen, und man 
erhält einen Fehler der näherungsweise proportional zur 
Nutzsignalleistung ist.

von Joe J. (j_955)


Lesenswert?

Detlef _a schrieb:
> Eine 2^x - sample FFT hat x 'stages', die jeweils die Multiplikation mit
> sin/cos machen und eine Addition.

Ich glaube pro Stage sind das N/2 kompl. MUL und eine kompl. ADD, wobei 
N die Anzahl meiner Punkte ist, die sich ergeben aus 2^x

Detlef _a schrieb:
> Bei 16Bit
> Eingangsdaten beträgt Dein SNR sowieso maximal ca. 96dB

Bedeutet das automatisch, dass ich mir jetzt das SNR ausrechnen kann und 
darüber eine Aussage habe, wie sich die Quantisierung auf mein 
Ergebnis(-> sprich Spektrum) auswirkt? Im Prinzip sind das ja dann nur 
die von Null versch. Werte, die im Spektrum auftauchen, richtig?

Meine mich zu erinnern, dass mir die Quantisierung beim AD/DA Wandler 
mal über den Weg gelaufen ist, kann ich dies dann soweit einfach 
übertragen und hier anwenden? Im Prinzip habe ich ja ein Signal und ein 
Quantisierungsfehlersignal...Hmm.

von Joe J. (j_955)


Lesenswert?

Jo-hann Ha-nnes schrieb:
> Im Prinzip sind das ja dann nur
> die von Null versch. Werte, die im Spektrum auftauchen, richtig?


Geimeint sind natürlich die Werte ausser den Spektralanteilen selbst....

von Detlef _. (detlef_a)


Lesenswert?

Zum Quantisierungsrauschen des AD-Wandlers:
http://de.wikipedia.org/wiki/Quantisierungsrauschen

Zum möglichen SNR bei integer-Zahlen: n-Bit erlauben Dir eine 
Quantisierung in 2^n Stufen, dh. Du machst einen 'Fehler' von 1/(2^n). 
Der 1 entspricht 0dB: 20*log10(1)=0. Der Fehler in dB beträgt 
20*log10(1/(2^n))= -n*20*log10(2) ~ -n*6.02 , also pi mal Daumen 6dB pro 
Bit.

Oder anders: Wenn Du einen Wert von 0 bis 1 mit einem Bit quantisierts, 
liegt Du im Mittel 0.5 daneben: 20*log10(0.5)=6.02dB SNR

Cheers
Detlef

von Joe J. (j_955)


Lesenswert?

Wie verhält es sich dann mit anderen Datentypen, glaube das mit dem 
Integer ist mir jetzt klar geworden, aber so kann ich das ja dann nun 
nicht auf Fließkommadatentypen übertragen, nicht?

von Theodor (Gast)


Lesenswert?

Da hat es ja doch nochgefunkt ;-)

Um mit Float oder Double alles durch zu rechnen reicht der Speicher 
meines MCs leider nicht aus.

Andreas Schwarz schrieb:
-----------------------------------------------------------------------
Da kannst du einfach float oder double verwenden und brauchst dir keine
Gedanken mehr über die Genauigkeit machen.

Ein Trick den man bei Festkommazahlen anwenden kann ist, die Daten vor
der FFT durch Schieben so zu skalieren dass der Wertebereich voll
ausgenutzt wird. Auf diese Weise kann man auch bei Signalen mit großem
Dynamikbereich oft mit relativ kurzen Wortlängen auskommen, und man
erhält einen Fehler der näherungsweise proportional zur
Nutzsignalleistung ist.
-----------------------------------------------------------------------

Kann ich das auch erreichen, indem ich vor den 
Multiplikationen/Additionen meine S16 auf double caste, die Rechnungen 
ausführe und dann wieder in S16 caste um den Wert abzulegen?
Gschwindigkeit spielt bei mir keine große Rolle.

Theodor

von Detlef _. (detlef_a)


Lesenswert?

Du weißt, dass float bzw. double beim AVR 4Byte lang sind?

Cheers
Detlef

von Theodor (Gast)


Lesenswert?

>Du weißt, dass float bzw. double beim AVR 4Byte lang sind?
Bei mir hat double 8 Byte und das beschränkt sich nicht auf den AVR.
Wie kommst du darauf, dass ich einen solchen verwende?

von Andreas S. (andreas) (Admin) Benutzerseite


Lesenswert?

Theodor schrieb:
> Kann ich das auch erreichen, indem ich vor den
> Multiplikationen/Additionen meine S16 auf double caste, die Rechnungen
> ausführe und dann wieder in S16 caste um den Wert abzulegen?

Klar, kannst du, wenn Laufzeit und Speicherverbrauch nicht stört.

von Andreas S. (andreas) (Admin) Benutzerseite


Lesenswert?

Jo-hann Ha-nnes schrieb:
> Wie verhält es sich dann mit anderen Datentypen, glaube das mit dem
> Integer ist mir jetzt klar geworden, aber so kann ich das ja dann nun
> nicht auf Fließkommadatentypen übertragen, nicht?

Fehler bei Floating-Point-Rechnungen lassen sich praktisch nicht 
berechnen, nur abschätzen. In der Regel geht man davon aus dass die 
Fehler bei double (64 Bit) so klein sind dass sie keine Rolle spielen. 
Um herauszufinden ob float (32 Bit) ausreichend ist führt man die 
gleiche Rechnung mit float und double durch und vergleicht das Ergebnis. 
Die Differenz ist das Quantisierungsrauschen das durch die Ungenauigkeit 
von float erzeugt wird, damit kann man das SNR für das Signal berechnen 
und entscheiden ob es ausreichend ist.

von Theodor (Gast)


Lesenswert?

>Klar, kannst du, wenn Laufzeit und Speicherverbrauch nicht stört.
In diesem Fall stört es nicht, ich habe ein Datenvolumen
von 16KB an 16-Bit-Werten im Speicher liegen, was dem halben RAM 
entspricht.
Ein paar hundert Byte zusätzlich für die Rechnungen mit double stören 
nicht,
ich kann nur nicht alles gleichzeitig mit double berechnen.
Danke, das hat mir sehr geholfen!

Theodor

von Joe J. (j_955)


Lesenswert?

@Autor:  Andreas Schwarz

thx!! Aber bei Integer ist es dann normales QUantisierungsrauschen, 
welches nur von meiner Bitbreite abhängt, wenn ich das richtig 
verstanden habe....Gut zu wissen!

von Helmut S. (helmuts)


Lesenswert?

> Ein paar hundert Byte zusätzlich für die Rechnungen mit double stören
> nicht,...

Wie soll das gehen?
In deiner FFT muss dann alles in double sein.
Du hast ein Riesenproblem mit deinem kleinen Speicher wenn du eine 
16k-FFT in long oder gar in double machen willst. Das mit int16 zu 
machen kann man vergessen.

von Theodor (Gast)


Lesenswert?

>Wie soll das gehen?
Indem nur der Butterfly mit double arbeitet.

>In deiner FFT muss dann alles in double sein.
Nein, warum den? Der Butterfly arbeitete vorher mit S16,
jetzt mit double, zweimal gecastet eben.

>Du hast ein Riesenproblem mit deinem kleinen Speicher wenn du eine 16k-FFT
>in long oder gar in double machen willst.
Es ist eine 8K-FFT mit 16-Bit-Werten.
Ich habe kein Problem damit, der Butterfly ist nun auf double geändert
und ich erhalte das richtige Spektrum zu meinen Testsignalen, wie
vorher auch, nur wird es jetzt wohl etwas genauer sein.

>Das mit int16 zu machen kann man vergessen.
Ich weiß nicht, was du für eine Implementierung verwendest, bei mir
klappt das super.

von Theodor (Gast)


Lesenswert?

P.S.: Wenn ich ein Signal mit 86db Dynamikbereich (16-Bit) reingebe, 
kann mein Spektrum auch höchstens diesem Bereich aufweisen.
Insofern macht die Verwendung von double für alle Ergebniswerte nicht 
viel Sinn. Da spätestens bei der 86-db-Schwelle der Rauschpegel liegt.

von Detlef _. (detlef_a)


Lesenswert?

Mein Beitrag vom 24.8.2010 ist fehlerhaft und ich muß ihn vollinhaltlich 
zurücknehmen. Die 6dB stimmen natürlich, aber die begründung war falsch. 
Richtig lautet sie so:


Ein Sinus 2*A Spitze/Spitze wird mit n Bit quantisiert, dann beträgt die 
Quantisierungsstufe q=2A/2^N . Der Quntisierungsfehler sei im Intervall 
-q/2 bis q/2 gleichverteilt, also ist die Wahrscheinlichkeitsdichte im 
Intervall konstant = 1/q .

Jetzt: Die Leistung eines stochastischen Signals der 
Wahscheinlichkeitsdichte p(e) beträgt :

sq=Integral -unendlich bis + unendlich von e^2*p(e)de, also bei uns:

integral -q/2 bis q/2 von e^2*1/q de =  (1/q) *(1/3)* q^3 von -q/2 bis 
q/2 =q^2/12

Der Effektivwert des Sinus beträgt A/sqrt(2), die mittlere Leistung also 
A^2/2

Mittlere Leistung im Verhältnis zu Rauschleistung also:

(A^2/2)/(q^2/12)=(A^2/2)/(4A^2/2^2N)=3/2*2^2N.

10*log(3/2*2^2N)= 1.76dB+N*6.02dB.

Schlüssel ist die Leistung eines stochastischen Signals der 
Amplitudendichteverteilung 1/q .

http://de.wikipedia.org/wiki/Quantisierungsrauschen
http://en.wikipedia.org/wiki/Quantization_error
https://ccrma.stanford.edu/~jos/mdft/Round_Off_Error_Variance.html

So ist's richtig
math rulez!

Cheers
Detlef

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.