Laden...

Klarheit betreffs FFT (Fast Fourier Transform) ?!

Erstellt von Nektarine vor 15 Jahren Letzter Beitrag vor 15 Jahren 3.991 Views
N
Nektarine Themenstarter:in
21 Beiträge seit 2008
vor 15 Jahren
Klarheit betreffs FFT (Fast Fourier Transform) ?!

Hallo liebe C#-Community!
Endlich habe ich mich auch angemeldet, nachdem ich schon länger passiver "Konsument" war. Zu meiner Schande muss ich gestehen, dass ich aber auch sofort eine Frage habe.

Wie bereits am Topic erkennbar - es geht um FFT.
Ich habe einige Sourcecodes mit verschiedenen FFT-Algos, aber so recht blicke ich da nicht durch, zumal auch scon alles optimiert ist und sich vom Grundalgorithmus stark entfernt hat. Wenn man mal bei Wikipedia FFT eingibt ist da ein Beispiel mit Pseudecode, das ich eigentlich noch recht durchsichtig finde, allerdings habe ich auch damit ein Problem. Da steht u.A. drin:

c(k) = g(k) + u(k) * e^(-2 * PI * i * k / n)

Soweit verstehe ich, dass es sich bei c(k) um das Array handelt, dass die Rekursion an die höhere Ebene zurück gibt bzw. g und u bereits von einer tieferen Rekursionsebene zurück gegeben wurden. Soweit, sogut. ABER:
was mache ich nur mit dem e^(...). Es handelt sich ja offensichtlich um komplexe Zahlen und u ist wohl der imaginäre Anteil. Aber wie kann ich das implementieren?! Ich würde es ja einfach so abtippen, aber das i ist ja die Wurzel aus (-1), also prinzipiel schwierig.
Vielleicht kann mir jemand talentierteres als ich es bin etwas klarheit verschaffen?!
Mich interessiert nur eine grundlegende FFT, nichts bis in die letzte Multiplikation optimiertes.
Und bevor ich es vergesse: Kann mir jemand erklären wozu bzw. wo ich bei dem Algorithmus ein Bytereverse-Methode brauche??!

Viele Fragen, ich weiß, sehr dreist für den ersten Beitrag, aber ich gelobe Besserung und hoffe, dass ihr mir trotzdem helfen könnt.

Vielen Dank und liebe Grüße,

Nektarine

Kaum macht man's richtig, schon funktioniert's!

1.361 Beiträge seit 2007
vor 15 Jahren

Hallo Nektarine,
willkommen als aktiver Part im Forum 😉

Es handelt sich ja offensichtlich um komplexe Zahlen

Du musst wirklich mit komplexen Zahlen rechnen.
Benötigst also einen neuen Datentyp (Klasse oder struct) "Complex" (der intern natürlich auch bloß aus zwei doubles/floats besteht, die ihrerseits Real- und Imaginärteil repräsentieren). Empfehlenswert ist wahrscheinlich eine fertige Implementierung zu benutzen, oder aber schnell selbst zu schreiben, musst ja nur die Addition, Multiplikation, ... auch noch als Funktionen einfügen.

Kann mir jemand erklären wozu bzw. wo ich bei dem Algorithmus ein Bytereverse-Methode brauche??!

Ich finde die Bezeichnung Bitreverse passender (schließlich wird ja nicht nur byte- sondern bitweise umgedreht)
Diese Invertierung ergibt sich durch die Aufteilung des Arrays in die geraden und ungeraden Stellen... immer und immer wieder.
Wenn man das bis zur untersten Ebene macht, hat man eine umsortierung des Ausgangsarrays, sodass die Indizes in Bitschreibweise genau umgekehrt sind.

Resultiert wohl daraus, dass das Zahlen im Dualsystem sind (2er) und bei jedem Schritt das Array auch in 2 hälften zerlegt, die darüberhinaus einmal am ende ne 0 haben (gerade) oder am ende ne 1 (ungerade).
Naja das ergibt sich eben gerade so 🙂
Schau dir mal die Grafik bei How the FFT works an.

beste Grüße
zommi

N
Nektarine Themenstarter:in
21 Beiträge seit 2008
vor 15 Jahren

Hallo! Danke erstmal für die schnelle Antwort. Die Sache mit dem Bit- und Bytereverse war ein Schreibfehler, meinte Bitreverse. Die Grafik auf der Seite konte mir hier auch einiges verdeutlichen.
Aber nochmal zu er Pseudocodezeile:

c(k) = g(k) + u(k) * e^(-2 * PI * i * k / n)

g ist also der Real- und u der Imaginärteil. Aber wenn ich ja nun für beide eine komplexe Zahl zurückbekommen habe, weiß ich immernoch nicht genau wie ich das implementieren soll. Wäre das dann einfach eine Addition/ Multiplikation?!
Nehmen wir an ich hatte eine Methode um 2 komplexe Zahlen zu Addieren / Multiplizieren, wäre das dann also einfach so:

c = complex.Addiere(g,u);  

bzw.

c = complex.Multipliziere(g,u);  

??

Ich bin offensichtlich kein großer Mathematiker =)

Kaum macht man's richtig, schon funktioniert's!

O
778 Beiträge seit 2007
vor 15 Jahren

Ja, das wäre es. (Auch wenn das Details sind, die mit der Mathematik nichts zu tun haben). Du kannst alternativ auch die Operatoren überladen, dann liest sich Code besser.

Ich kenne FFT nicht, aber warum geht ihr alle davon aus, dass g real und u imaginär ist? e^(-2Pii*k/n) ist nichts weiter als ein Vorzeichen. Das kann aber abhängig von k und n alles auf dem Einheitskreis sein. Bei k=0 oder k=n wäre das Vorzeichen z.B. ein + und zwar ein stinknormales reales +.

N
Nektarine Themenstarter:in
21 Beiträge seit 2008
vor 15 Jahren

Danke für deine Antwort.

Ja, das wäre es. (Auch wenn das Details sind, die mit der Mathematik nichts zu tun haben). Du kannst alternativ auch die Operatoren überladen, dann liest sich Code besser.

Dass das mit dem Quelltext so wäre weiß ich, meine Frag zielte aber doch auf das technische Detail ab, und zwar ob es denn nun z.B. eine Addition ODER eine Multiplikation wäre. Für mich ist das ein etwas böhmischer Wald. Und die Sache mit dem Vorzeichen stimmt so glaube ich (mit meinem mathematischen Verständnis) nicht:
Der Term (-2Pii*k/n) ist schließlich ein Exponent und kann daher sehr starken Einfluss auf das Gesamtergebnis haben.
Einfaches Beispiel: Ist k = 0, dann ist e^(...) = 1; (das komplexe "i" mal außer Acht gelassen 🤔 )
Ist k = n, dann ist e^(...) = 0.0019; (wieder ohne i)
Der Unterschied ist also beträchtlich, erstrecht bei großen u(k).

Vielleicht ist die Sache mit der Pseudocodezeile auch ein falschen Ansatz.
Unter dem Link von zommi (http://www.dspguide.com/ch12/2.htm) habe ich gelesen, was eigentlich offensichtlich ist: c(k), u(k) und auch g(k) sind komplexe Zahlen, es geht also nur darum u und g zu c zu kombinieren. Hierzu muss das irgendwie mit einer Sinusschwingung (=>"sinusoid") multipliziert werden, was glaube ich der e^(...)-Teil ist, allerdings habe ich keine Ahnung wie. Aus dem Fortran-Code eine Seite weiter (im Dokument von zommis Link) kann ich es irgendwie nicht entnehmen...

Kaum macht man's richtig, schon funktioniert's!

J
1.114 Beiträge seit 2007
vor 15 Jahren

Vorsicht: e^(ix) ist was gaaaaaaanz anderes wie e^x.

Der imaginäre Exponent mach nun nämlich was ganz Spezielles aus der Exponentialfunktion, nämlich:

e^(ix) = cos(x) + i*sin(x)

Das kannst du für deine FFT nutzen. e^(ix) ist nur eine andere Schreibweise.

Aber um die grundlegenden komplexen Rechenoperation (Addition und Multiplikation) kommst du nicht rum.

1.361 Beiträge seit 2007
vor 15 Jahren

Hi Nektarine,

das mit der Fourier-Transformation ohne großes Hintergrundwissen ist schwierig 😉
Hier ist übrigens noch ein Link zur Uni Flensburg. Besonders dann der Part mit dem Quellcode könnte wieder interessant sein. wobei es dort natürlich auch schon wieder optimiert wurde.

Vielleicht solltest du, vor allem fürs Verständnis, erst anfangen ne stupide Fouriertransformation zu implementieren. Keine FFT.
Heißt: keine Rekursion, keine halbierung, kein O(nlogn) Aufwand, sondern O(n^2).

Was musst du für diese Basic-Variante nun haben?1.ein wenig Kenntnis von komplexen Zahlen und deren Addition/Multiplikation. (schau dir insbesondere die Darstellung in der Gaußschen Ebene an und in dem Zusammenhang den Einheitskreis mit den nten Einheitswurzeln (unityRoot)) 1.eine Implementierung von komplexen Zahlen in C#
Hierzu bietet sogar die MSDN ein schönes Beispiel mit Operatorüberladung.
Lediglich die Multiplikation fehlt, aber mit Hilfe von Wiki solltest du die ohne Probleme hinzufügen können.

1.eine Funktion, die dir e^(ix) berechnet. Denn wie schon Jelly erwähnt hat, passiert bei komplexen Argumenten (und i ist komplex) etwas ganz anderes als im rein reellen. Gleichzeitig hat Jelly schon eine andere Darstellung angegeben, die es dir (mit Hilfe von Sin und Cos aus System.Math) ermöglicht leicht das komplexe Ergebnis zu berechnen. Und darauf aufbauend dann eine Funktion

expi(double k, double n)=exp(-2*Pi*i*k/n)=new Complex(cos(-2 * Pi * k/n) , sin(-2 * Pi * k/n))

zu erstellen. 1.Die Grundformel für komplexe DFT

Letztere Summe lässt sich somit für dein Array mit zwei verschachtelten Schleifen berechnen (daher dann auch O(n^2)):


Complex[] DFT(Complex[] a)
{
 int n = a.Length;
 Complex[] b = new Complex[n];
 for(int k=0; k<n; k++)
 {
  b[k] = new Complex(0,0);
  for(int j=0; j<n; j++)
  {
  b[k] += expi(j*k, n) * a[j];
  }
 }
 return b;
}

Sollte hoffentlich funktionieren, wenn ich mich jetzt nicht vertan hab. 😉
Soweit solltest du erstmal genug zum experimentieren haben 🙂
Wenn das dann wirklich läuft, dann fang an mit der FFT.

beste Grüße
zommi

1.361 Beiträge seit 2007
vor 15 Jahren

Desweiteren wird bei

c(k) = g(k) + u(k) * e^(-2 * PI * i * k / n)

mit g(k) und u(k) nicht Real- und Imaginärteil bezeichnet, sondern es ist jeweils das k-te ungerade Element, bzw. das k-te gerade Element des Arrays gemeint, die ihrerseits beide komplex sind.
Sprich {a[0], a[1], a[2], a[3], ...} = {g[0], u[0], g[1], u[1], ...}.

//Edit:
Sorry, wollte editieren und nicht anfügen

O
778 Beiträge seit 2007
vor 15 Jahren

Abhängig davon, was du mit den komplexen Zahlen machst kann es auch sinnvoll sein nicht die Koordinaten- sondern die Polarform der komplexen Zahlen zu speichern. In der Polarform lassen sich komplexe Zahlen deutlich einfacher multiplizieren und e^(ix) verkommt zum Witz (r=1 und phi=x). Einzig kompliziert ist die Addition (und damit logischerweise auch die Subtraktion). Aber das sind alles wieder Optimierungsideen...

N
Nektarine Themenstarter:in
21 Beiträge seit 2008
vor 15 Jahren

Ich muss mich für die tatkräftige Unterstützung bedanken. Kurz bevor ich zommi's Post las, hab' ich dann doch nochne FFT hinbekommen, die mir aber irgendwie eigenartige Werte brachte. Dann hab ich also die DFT wie Zommi sie gepostet hat probiert und bekam ganz ähnliche Werte, alsodenke ich mal, dass sie wohl doch stimmen. Aber die Werte sind trotzdem recht eigenartig:
Sie wandern von x bis -x zurück zu x bei einem Feld der Länge 16, wo ich abwechselnd als real-Teil x und 0 eingetragen habe.
Das Ergebnis ist völlig symmetrisch. Bekomme ich also nur halb soviele Frequenzbänder raus?

Übrigens toll, der letzte Kommentar auf der MSDN-Seite:

Java unterstützt das Überladen von Operatoren nicht, obwohl es intern den +-Operator für Zeichenfolgenverkettungen überlädt.

...immer diese Javaprogrammierer.... 😜 👅

Kaum macht man's richtig, schon funktioniert's!

1.361 Beiträge seit 2007
vor 15 Jahren

Hallo Nektarine,

poste doch mal dein Eingabe-Array. So ganz klar ist mir das aus deiner Beschreibung nicht 🤔

Aber eine gewissen Symmetrie wird es immer bei reellen Eingangssignalen geben.
Wie sonst werden aus n rellen Zahlen als Ergebnis n komplexe Zahlen (die ja intern quasi aus zwei reellen Zahlen bestehen). Wo soll sich die Information verdoppeln?

Das ergebnis eines reellen Inputs ist hermitesch.
Das heißt, dass die realteile symmetrisch sind und die imaginärteile antisymmetrisch.
also re(-k)=re(k) und im(-k) = -im(k).

Demzufolge sind die n komplexen Zahlen untereinander abhängig und es gibt wieder nur n Freiheitsgrade (alles wieder gerettet 😁 )

Vielleicht hast du fläschlicherweise bei deiner Symmetriebetrachtung auch nur die Realteile angeschaut? Oder aber auch die Imaginärteile sind symmetrisch; dann müssen sie aber notwendigerweise alle 0 sein.
Dies würde bedeuten, es gibt nur Cosinus-Anteile im Eingangssignal aber keine Sinusanteile. (hoffentlich hab ichs jetzt nich wieder vertauscht 😉)

PS:

und bekam ganz ähnliche Werte

ganz ähnliche != identische
Das sollte mindestens einem von uns beiden zu denken geben. 😉

N
Nektarine Themenstarter:in
21 Beiträge seit 2008
vor 15 Jahren

Hallo!

Als Engangsarray hatte ich folgendes:


complex[] arr = new complex[16];
for (int i=0; i<8; i++)
{
    arr[2*i]= new complex(1,0);
    arr[2*i+1] = new complex(0,0);
}

Und zum ganz ähnlich:Bei der FFThatte ich noch einen normierungsfaktor drin, der das auf Werte zwischen -1 und 1 genormt hat, waren also gleiche Werte, nur normiert.

Ich glaube ich habe die Sache mit dem imaginärteil nicht richtig gemacht, und zwar habe ich als Frequenzwert einfach den Betrag der komplexen Zahl genommen, daher wahrscheinlich die Symmetrie. Aber was ist denn nun eigentlich der echte Frequenzwert?!
Ganz schön schwierig, aber ich bin schon ganz schön weit gekommen - durch eure Hilfe! 🙂 🙂 🙂

Ach und noch eine Frage: Wie kann ich eigentlich die dB-Werte ausrechnen? Muss das Programm quasi wissen, dass es sich beispielseise um eine 16Bit Wavedatei handelt und kann dann bei voller Aussteuerung sagen 96dB und bei niedrigeren Werten entsprechend runter rechnen oder gibts da eine bessere Möglichkeit?!

//Edit: Tippfehler

Kaum macht man's richtig, schon funktioniert's!

N
Nektarine Themenstarter:in
21 Beiträge seit 2008
vor 15 Jahren

Hab meine letzte Frage betreffs der dB selbst lösen können:

der dB-Wert eines Samples ist 20*log(Samplewert / 2^(Auflösungin Bit))

//Edit: bzw. ist das der Abstand in dB von der maximalen Aussteuerung.
Maximale Aussteuerung = 20*log(2^(Auflösung in Bit))

Kaum macht man's richtig, schon funktioniert's!