FFT ohne Floatingpoint

[Charly_cs]

Hallo!

Hab mal wieder ein wenig Zeit mich mit der FFT zu beschäftigen und wollte diese für den mega644 in C programmieren jedoch wie der Titel schon sagt ohne Floatingpoint.
Möchte folgendermaßen vorgehen:
Da ich ja weiß das ich 1024 Werte aufnehme kann ich die benötigten Drehfaktoren schon vorher berechnen und in einer Tabelle ablegen. Entweder programmiere ich sie mir als Konstante oder lasse diese Tabelle bei einschalten des µC sozusagen als Initialisierung vorrechnen. Die eigentliche FFT-Funktion pickt sich dann die benötigten Ergebnisse raus usw.
Somit gibt es keinen Sin/Cos Berechnungen mehr, sondern Addition, Multiplikation und Subtraktion. Da das alles allerdings mit Floatinpoint laufen würde, wäre der µC seeehr langsam... Bräuchte jetzt einen Tip wie ich das ganze mit Ganzzahltypen also Integern lösen kann.
Wenn ich die Drehfaktorentabelle eben mit integern füllen würde, wäre der restliche Weg klar... Meine Idee ist z.B. den
Sin(0.5) = 0,479425538 rad
mit 1000 zu multiplizieren und den Rest hinter dem Komma abzuschneiden: 0,479425538 * 1000 = 479 --> Ergebnis ist ein Integer.
So wollte ich mit allen anderen Sin/Cos Berechnungen vorgehen. Am Ende wäre die Tabelle nur noch mit Integern gefüllt. Denke ich da in die richtige Richtung oder ist das falsch?

Grüße
Charly

[mare_crisium]

Charly_cs,

ja, Dein Weg ist gangbar. Du musst nur aufpassen, dass die Integers lang genug sind, um die grösste vorkommende Zahl ohne Überlauf aufnehmen zu können.

Es gibt noch einen anderen Weg: Nimm einfach die Fast-Hadamard-Transformation (FHT) und rechne das Ergebnis anschliessend in Fourier-Koeffizienten um. Das ist vorteilhaft, weil die FHT nur mit Subtraktionen und Additionen und ohne Drehfaktor-Tabelle auskommt. Weil die Hadamard-Funktionen (sind eigentlich Walsh-Funktionen) ein orthogonales Funktionensystem bilden, ist die Abbildung auf den Fourierkoeffizienten ein-eindeutig. Du kannst also ohne Informationsverlust hin- und zurücktransformieren. Deshalb ist auch die Transformationsmatrix Hadamard->Fourierkoeffizienten symmetrisch. Es genügt also, die obere rechte Dreiecksmatrix zu berechnen. Ich würde die Werte im Voraus auf dem PC berechnen und dann als Tabelle im Programmspeicher des AVR ablegen.

mare_crisium

[Charly_cs]

Hallo!

Danke für die schnelle Antwort! Auf das Problem mit dem Überlauf bin ich auch schon gekommen. Werde da einfach durch rumprobieren die max. Werte finden und somit den Überlauf ausschließen. Simuliere den Code in LabView.
Zum FHT: Habe davon leider noch nie etwas gehört. Kannst du mir bitte eine Quelle nennen, wo ich genaueres darüber nachlesen kann?
Danke!

Charly

[mare_crisium]

Charly,

das mit dem Überlauf würde ich nicht mit Trial und Error ermitteln. Ich, jedenfalls, habe die starke Tendenz, die Randbedingungen und Einschränkungen meiner Algorithmen schnell zu vergessen. Und wenn's dann doch mal passiert, und eine Einschränkung überschritten wird (in Deinem Fall: Ein Datenpunkt den Überlauf verursacht), dann such' ich mir einen Wolf, bevor ich mich daran erinnere, dass ja der Überlauf auftreten könnte. - Also rate ich Dir, lieber ein Stück Papier zu nehmen und zu versuchen, bei vorgegebener Integerlänge abzuschätzen, wie gross der maximal zulässige Betrag eines Datenwertes sein darf. Dann kannst Du eine Fehlermeldung erzeugen, wenn die Einschränkung verletzt ist, die Dir sehr viel Suchzeit erspart ;-).

Leider ist die Hadamard-Transformation noch nicht so populär. Hauptsächlich wird sie in der Akustik, der Übertragungstechnik (zellulare Funknetze) und der Kyptographie verwendet. Ich habe aber auch schon Anwendungen aus der Regelungstechnik gesehen. Für die Bildverarbeitung könnte sie auch nützlich sein.

Hier ist eine kleine Linksammlung:
http://www.module.ru/files/nm6403fht-e.pdf
Beschreibung der Hadamard-Transformation und der FHT Seite 4ff

http://www.cs.technion.ac.il/Labs/Isl/Project/Projects_done/VisionClasses/DIP/Hadamard_Walsh/node3.html
Kurze Einführung in die Hadamard-Walsh-Transformation

http://www.cs.technion.ac.il/Labs/Isl/Project/Projects_done/VisionClasses/DIP/Hadamard_Walsh/node5.html
C-code zur FHT

http://www.akustik.rwth-aachen.de:9673/pub/material/messtechnik.pdf
Sehr kurze Einführung der FHT in Zusammenhang mit Raumschall und MLS Seite 55ff

http://www.udobarth.de/Profil/Doku/Uni/Diplom.pdf
Kurze Einführung in die Hadamard-Walsh-Funktionen und die FHT Seite 16ff

mare_crisium

[Charly_cs]

hallo!

Vielen Dank für die Linksammlung! Habe mir alles angeschaut und bin ein wenig verwirrt. Wie rechne ich die Ergebnisse der FHT in Ergebnisse der FFT um?
In der FHT wir ja der Drehfaktor immer 1 bzw. -1 angenommen wodurch die Multiplikation entfällt und damit auch die spektralen Anteile des Sin/Cos. Damit kommt bei der FHT ein komplett anderes Ergebnis raus. Hab das ganze mal mit einer 4Punkte FFT und FHT gemacht (Anhang).

Gruß
Charly

[robocat]

ist mir bissl peinlich das zuzugeben, aber ich verstehe diese fft immer noch nicht so ganz. dein rechenbeispiel ist schonmal sehr aufschlussreich, aber was sind die eingangswerte (1,2,3,0 oder?) und was sagt das ergebnis aus?

@topic: ich habe für ein projekt die für avr-µC optimierte fft von elmar chan (hier zu finden: http://elm-chan.org/works/akilcd/report_e.html) übernommen. ich glaube dort werden fixedpoint-zahlen verwendet.

gruesse und danke wenn sich jemand die mühe macht, kurz ein wenig zu erklären.

[mare_crisium]

@Charly,

ich schreib' dazu was auf - dauert aber noch ein bisschen. Im Prinzip bekommst die Transformationsmatrix von Hadamard- zu Fourierkoeffizienten, indem Du die beteiligten Winkelfunktionen Hadamard-transformierst.

@robocat,

das braucht Dir überhaupt nicht peinlich sein! Is' ja auch nicht das einfachste Thema.

Bezieht sich Deine Frage auf das Spreadsheet von Charly? In dem Beispiel hat er ausprobiert, wie sich die Ergebnisse von Fourier- und Hadamard-Transformation voneinander unterscheiden. Die Eingangswerte hat er willkürlich ausgewählt. Natürlich kommen dabei zwei ganz verschiedene Sätze von Zahlenwerten heraus. Man kann die Ergebnisse aber zwischen den beiden Systemen aber hin- und her "übersetzen". Oder fragst Du wegen der Grundsätze der Analyse nach orthogonalen Funktionensystemen? Vllt hilft dann der Anhang.

mare_crisium

Edit: Attachment gelöscht wegen Upload-Quota. Bei Interesse pn schicken!

[robocat]

danke mare_crisium, das pdf ist sehr verständlich geschrieben.

ich hatte gehofft, dass man das rechenbeispiel von charly für die fft in einem real-world beispiel erklären kann. zB angenommen es wären 4 samples, je 1 pro 1000stel sekunde, dann müsste ich doch über die fft 2 frequenzkomponenten bekommen, den anteil von 0-500Hz und den von 500-1000 (oder?). aber er berechnet 4 werte, und in welchem zusammenhang zB der wert 6 zur "amplitude" der zugrundeliegenden werte steht, sehe ich auch nicht. oder verstehe ich etwas völlig falsch?

ich denke aber, ich sollte mir erst nochmal die zugrundeliegende theorie ansehen, bevor ich das ganze in einem praktischen beispiel erklärt haben möchte.

danke für deine antwort & gruesse von der katze

[mare_crisium]

So, Charly,

hier ist die Geschichte mit der Umrechnung von Hadamard- in Fourierkoeffizienten. - Wie üblich hab ich's eher populär- als exakt-mathematisch aufgeschrieben, in der Hoffnung, dass die Sache so auch für eventuelle Mitleser verständlicher wird.

@ robocat,
Wenn Du also 4 Messwerte hast, bekommst Du auch 4 Fourier-Koeffizienten. 2 davon stellen die Amplitude der Schwingungen dar, die anderen beiden ihre Phase. Wie gross der Beitrag der Schwingungen ist, deren Periode nicht mit den Grundfrequenzen übereinstimmt, überblicke ich nicht. Das müsstest Du selbst mal ausprobieren.

mare_crisium

Edit: Attachment gelöscht wegen Upload-Quota. Bei Interesse pn schicken!