FFT: Numerische Genauigkeit?

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

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

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

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.

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.

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

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

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?

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

Cheers
Detlef

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.

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.

@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!

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

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