Blocktransformation beliebiger Größe

PhilippHToner

Hallo miteinander,

ich stehe einmal wieder vor einem Rätsel. Ich baue ein Effekt-Plugin für eine Musiksoftware. Die Musiksoftware liefert N Samples und ich würde gerne das Frequenzspektrum anzeigen lassen. Ich kann leider nicht davon ausgehen, dass ich den cooley-tukey fft hernehmen kann (den ich schon programmiert habe), da nicht definiert ist, wie viele Samples eingehen.

Nach kurzer Recherche würde der Bluestein-FFT funktionieren. Durch geschickte Umwandlung des DFT entsteht eine Faltung, die man wiederum anhand des Faltungstheorems umwandeln kann wie es auf Wikipedia beschrieben steht:
[url]http://en.wikipedia.org/wiki/Bluestein's_FFT_algorithm[/url]
http://de.wikipedia.org/wiki/Faltung_(Mathematik)#Faltungstheorem

"[*]" = math. Faltung
"F(...)" = math. fast fourier trans.
"iF(...)" = math. inverse fast fourier trans.
"e" = math. Element von

Also entsteht folgende Faltung: X_k = b_k^* * (a_n [*] b_n)

Durch das Faltungstheorem entsteht dann:

M=2^r, r e Z, M >= N
A_n = { a_n, 0, ..., 0 } Länge M
B_n = { b_n, 0, ..., 0 } Länge M

X_k = b_k^* * iF( (2π)^M/2 * F(A_*) * F(B_*) )_k

Fragen/Probleme:
- ist das so richtig?
- (2π)^M/2 wird für z.B. 1023 nicht mehr kalkulierbar

vielen Dank
gruß Philipp

lustig

Ich verstehe dein Problem noch nicht ganz, könntest du das vielleicht nochmal genauer beschreiben?

PhilippHToner schrieb:

Die Musiksoftware liefert N Samples und ich würde gerne das Frequenzspektrum anzeigen lassen. [...], da nicht definiert ist, wie viele Samples eingehen.

Wieso ist das ein Problem? In welcher Grössenordnung liegt denn N? Und was heisst "nicht definiert"?
Ist das Problem, dass du N nicht als 2^k schreiben kannst, für Cooley–Tukey? Was spricht denn gegen zero padding?

Ausserdem du scheinst eine mehr oder weniger zufällige Formel aus dem Wiki Artikel zur Fouriertransform zu verwenden. Wieso glaubst du, dass die auch für die DFT gilt? (Lies z.B. mal http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Circular_convolution_theorem_and_cross-correlation_theorem). Die diskrete Faltung kann man durchaus mit der FFT effizient implementieren ( http://en.wikipedia.org/wiki/Convolution#Fast_convolution_algorithms)... aber dazu solltest du dich vielleicht mal mit den Eigenschaften der DFT bezüglich Faltung / zyklischer Faltung auseinandersetzen.

PhilippHToner

Okay, sorry dann nochmal genauer:

Ich habe N reelle Eingangswerte (Samples aus einem Lied z.B.). Das sind einfach die chronologischen diskreten Ortswerte zwischen [-1;1]. Größe von N ist unbekannt aber wahrscheinlich unter 4000 Werten. Diese Werte bekomme ich immer häppchenweise. Jedes Mal, wenn ich wieder welche bekomme, möchte ich das Frequenzspektrum anzeigen lassen, bis wieder neue kommen. Wie z.B. ein Mediaplayer seine Visualisierung füttert, die wiederum auf die Musik reagiert, z.B. bei tiefen Frequenzen zuckt oder Farbe wechselt.

Kurz gesagt nochmal: ich brauche eine BLocktransformation mit beliebigem N (wird nicht zu groß)

Problem:
- cooley-tukey geht nur für 2^k Elementen (zero padding verfälscht das Ergebnis doch; angenommen ich bekomme 2049 Elemente, muss ich bis 4096 padden)
- bluestein kann ich nachvollziehen, aber es kommen zu hohe Werte vor, mit der der Computer nicht mehr rechnen kann

Ich bräuchte wie der Titel schon sagt, eine Blocktransformation beliebiger Länge.

Oder Folgende Überlegung:
Bestünde die Möglichkeit anstatt mit Nullen, einfach mit den Samples von vorne wieder zu padden, bis ich zu 2^k Elementen komme?

Ich habe z.B. N=10
Samples abcdefghij => abcdefghijabcdef (N=16)

gruß Philipp

krümelkacker

Wenn Du eine schnelle Faltung per FFT implementieren willst, schaust Du Dir am besten von diesem Buch http://dspguide.com/ mindestens Kapitel 18 an. Die Kapitel kann man sich da sogar als PDF frei runterladen.

Der Unterschied zwischen einer kontinuierlichen Fourier-Transformation und einer diskreten ist der, dass Du mit der diskreten nur eine "zyklische Faltung" berechnen kannst, was zu einem "Wrap-Around"-Fehler führt, wenn man das einfach ignoriert.

Einen Hall-Effekt habe ich per "FFT-Faltung" auch schon gebaut. Ich weiß also, wovon ich rede.

PhilippHToner

Ja bei mir wird der Faltungshall die nächste Baustelle werden. Ja ich Idiot vergesse immer das Buch. Ich hab ungefähr schon auf 100 Fragen das Buch als Antwort bekommen. Ich werde mir des Ding komplett reinpressen, damit ich hier nicht immer soo viele Fragen stellen muss und auch mal beantworten kann.

Merci, ich schreibe, sobald ich die Lösung weiß/vermute

gruß Philipp

PhilippHToner

Ich bin wieder zurück :-D. Ich habe alles gelesen zzgl. einem gekauften Buch. Ja also wir springen gleich von der Transformation zum Faltungshall. Ich habe meine Eingangswerte x und meinen Faltungskernel h. Bei der Lösung von x*h als Faltung bin ich über die Möglichkeit der Overlap-Add, bzw. Overlap-Save gestolpert. Diese sind bestimmt sehr toll für "kleinere" Kernel. Mit klein meine ich nicht in der Zahl der Samples, sondern in der zeitlichen Länge, die diese darstellen. Was ist, wenn ich mit 20sec Hall arbeiten möchte. Dann müsste bei 44,1 Samples/sec doch 20*44100 Sample Kernels entstehen.

Es gibt ja die Lösung mit einer Tapline Länge L, die ein h simuliert, dass an erster Stelle und nach L Samples einen Impuls haben. Aber mehr schon auch nicht mehr.

Wenn ich mein h fenstere, dann geht mir doch die Info verloren, dass sich das Frequenzspektrum von h zeitlich ändert. Liege hier richtig oder habe ich einen Denkfehler?

In Büchern etc. reden die immer über h mit M knapp größer als 60 Samples wo sich dann FFT-Convolution erst wirklich rentiert. Aber Ordnung 40000 Samples nie.

Vielen Dank für die Antworten

krümelkacker

PhilippHToner schrieb:

Wenn ich mein h fenstere, dann geht mir doch die Info verloren, dass sich das Frequenzspektrum von h zeitlich ändert. Liege hier richtig oder habe ich einen Denkfehler?

Das ergibt für mich so keinen Sinn.

PhilippHToner schrieb:

In Büchern etc. reden die immer über h mit M knapp größer als 60 Samples wo sich dann FFT-Convolution erst wirklich rentiert. Aber Ordnung 40000 Samples nie.

Wenn es sich ab 60 lohnt, dann ganz besonders bei 40000, nicht? Zumindest in der Theorie.

PhilippHToner

Hmm irgendwas stimmt noch nicht.
Ein Problem bei Hall ist doch, dass Nachmodellierungseffekte, wie sie in alle gängigen Audioprogrammen zu finden sind, sehr primitiv arbeiten und meist einen Runden Raum simulieren mit Parameter wie Dämpfung (Hi/Lo), erste Reflektionen, Länge. Klar gibt es komplexere, aber wenn man einen Raum 1:1 als Verhallungsgrundlage benutzen will (und ihn nicht 3D nachmodelliert und mit Raytracing simuliert :D) braucht man die Hallantwort des Raumes. Und die bekommt man doch, wenn in dem Raum ein Weißes Rauschen abgespielt und aufgenommen wird.

Wenn ich nun möchte, dass ein Gesang in diesem Raum simuliert wird, falte ich den Gesang mit der Hallantwort h aus diesem Raum. Ist doch nichts anderes wie ein FIR. Theoretisch könnte man h auch mit Polen und Nullstellen approximieren.
Nachdem die Hallantwort ziemlich lang sein kann (z.B. Lagerhalle), kann mir passieren, dass ich neben meinem (unendlich langen) x auch ein sehr langes h bekomme. Und diese Overlap-Methoden basieren ja darauf, dass bei einer Faltung einer der beiden Operanden (x oder h) zu lang bzw. der anderen kurz genug ist, um eine Faltung im Spektralbereich durchzuführen.
Deshalb die Frage, gibt es eine Overlapmethode, die sowohl lange x, als auch lange h mit kleineren ffts darstellt?

Beste Grüße

krümelkacker

gibt es eine Overlapmethode, die sowohl lange x, als auch lange h mit kleineren ffts darstellt?

Eine Faltung ist linear und verschiebungsinvariant. Von daher kannst Du h mit
h = [h1 h2 h3 ... hn] in n Blocke zerlegen, x mit separat mit h1,h2,... falten und die Einzelergebnisse hinterher per Addition zusammen setzen. Du reduzierst mit kleinerer FFT-Blocklänge die Input/Output-Latenz (angenommen, man rechnet unendlich schnell), aber der Rechenaufwand wird auch größer. Trotzdem kann es so auch schneller klappen, weil die Speicherzugriffe da anders aussehen und eventuell Cache-freundlicher sind. Ich würde also die größtmögliche FFT-Länge nehmen, die noch vertretbar ist (bzgl Cache und Latenz) und den Rest "zu Fuß" ausrechnen.

Außerdem fällt mir da gerade noch ein Trick ein, mit dem Du die Sache noch etwas beschleunigen kannst. Die einzelnen "Spektralblöcke" FFT(h_k) kannst Du vorberechnen. Anschließend kannst Du dann zwei Faltungen mit nur einer FFT und einer iFFT berechnen:

H1 = FFT(h1);
...
x1 = ...
x2 = ...
zz = iFFT( FFT(x1 + x2 * 1i) .* H1 );
f1 = real(zz); // = conv(x1,h1)
f2 = imag(zz); // = conv(x2,h1)

sofern x und h beide reellwertig waren.

PhilippHToner

Hi,

danke für die Info, gut dass es geht :).

Das mit dem Vorberechnen hätte ich auch schon im Hinterkopf damit ich mir H# = FFT(h#) sparen kann.

Ich implementiere mit OpenCL und denke ich werde die Overlap-Save Methode verwenden, damit ich ohne zu syncen zurückschreiben kann. Hier muss ich allerdings nochmal recherchieren, wie ich welche Blöcke verrechnen muss. Angenommen meine h#-Blöcke haben Länge L# und meine x#-Blöcke haben Länge M#.

x# und h# sei ein Block zu berechnen mit y# = iFFT( FFT(x#) * FFT(h#) ).
Dann muss x# und h# die gleiche Länge haben. Ich berechne so doch ein y#, wobei ich die "rechten" Samples nicht brauche. Die Frage die ich noch habe ist, wie groß ist "rechts" und wie groß ich M# bzw. L# wählen muss, dass Overlap-Save funktioniert. Ich möchte schließlich FFT mit tukey-cooley Big-O(N*log(N)) berechnen und muss sowieso zeropadden, wie es Overlap-Save ohnehin macht.

Wikipedia habe ich mir schon angeschaut, aber mit den h# Blöcken komme ich jetzt durcheinander.

Vielen Dank!