😈 🤴🏻 👊🏻 Fourier-Transformation. Das schnelle und das wütende 🎆 🤛🏾 🛌🏾

Bei der Entwicklung von Algorithmen stoßen wir häufig an die Grenze der Rechenkomplexität, die anscheinend nicht zu überwinden ist. Die Fourier-Transformation ist komplex

O (n^{2})

$O (n ^ 2)$ und eine schnelle Version, vorgeschlagen um 1805 von Haus ¹ (und 1965 von James Cooley und John Tukey neu erfunden)

O (n l o g (n))

$O (nlog (n))$ . In diesem Artikel möchte ich Ihnen zeigen, dass Sie die Konvertierungsergebnisse in linearer Zeit erhalten können

O (n)

$O (n)$ oder sogar konstante Schwierigkeit erreichen

O (1)

$O (1)$ unter bestimmten Bedingungen, die in realen Problemen gefunden werden.
Fourier-Analyse

Als ich vor der Aufgabe stand, ein Programm zur Analyse der Übertragungsfunktionen von Soundsystemen in Echtzeit zu schreiben, wandte ich mich wie alle anderen zunächst der schnellen Konvertierung zu. Alles war in Ordnung, aber mit der Größe des Zeitfensters wurde die CPU-Auslastung unangemessen groß und es musste etwas getan werden. Es wurde beschlossen, die Transformation anzuhalten und erneut zu studieren und gleichzeitig nach Wegen zu suchen, um das Problem zu lösen. Kehren wir zur ursprünglichen Transformation von Joseph Fourier ^{2 zurück} :

f (x) = s u m l i m i t s_{- i n f t y}^{+ i n f t y} c_{k} e^{2 p i i k x / T} c_{k} = f r a c 1 T i n t l i m i t_{0}^{T} f (x) e^{- 2 p i i k x / T} d x

$f (x) = \ sum \ limits _ {- \ infty} ^ {+ \ infty} c_ke ^ {2 \ pi ikx / T} \\ c_k = \ frac {1} {T} \ int \ limit_0 ^ Tf ( x) e ^ {- 2 \ pi ikx / T} dx$

Wir werden uns genau ansehen, was hier los ist. Jeder Ausgabewert im Frequenzbereich

c_{k}

$c_k$ ist die Summe aller Eingangswerte des Signals

f (x)

$f (x)$ multipliziert mit

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ . Um Berechnungen durchzuführen, müssen wir alle Eingabedaten für jeden Ausgabewert durchgehen, d. H. diese zu erfüllen

n^{2}

$n ^ 2$ Operationen.

N loswerden

Ich möchte Sie daran erinnern, dass es ursprünglich darum ging, Schalldaten in Echtzeit zu analysieren. Dazu wird das ausgewählte Zeitfenster (im wesentlichen ein Puffer) der Größe N mit Daten mit einer Frequenz f _d gefüllt, die der Abtastfrequenz entspricht. Mit der Periode T werden die Eingabedaten vom Zeitfenster in das Frequenzfenster konvertiert. Wenn Sie sich die reellen Zahlen ansehen, variiert N von 2 ¹⁴ (16 384) bis 2 ¹⁶ (65 536) Stichproben (die Werte werden von der FFT geerbt, wobei die Fenstergröße eine Zweierpotenz sein muss). Zeit T = 80 ms (12,5 Mal pro Sekunde), wodurch Sie Änderungen ganz bequem sehen und die CPU und die GPU nicht überlasten können. Die Abtastfrequenz f _{d ist} Standard und beträgt 48 kHz. Berechnen wir, wie viele Daten im Zeitfenster sich zwischen den Dimensionen ändern. Während der Zeit T tritt es in den Puffer ein

48000 * 0, 08 = $ 384

$48000 * 0,08 = $ 384$ Proben. Somit werden nur 5% bis 23% der Daten im Fenster aktualisiert. Im schlimmsten Fall fallen 95% (und bestenfalls 73%, was auch viel ist!) Der verarbeiteten Proben immer wieder in die Konvertierung, obwohl sie bereits in früheren Iterationen verarbeitet wurden.

Der aufmerksame Leser in diesem Moment wird seine Hand heben und sagen: „Warte, aber was ist mit dem Koeffizienten?

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ ? Schließlich werden bei jeder neuen Transformation dieselben Daten an den neuen Positionen der Reihe lokalisiert und haben daher unterschiedliche Koeffizienten? “ Erinnern wir uns für alle fünf für ihre Pflege an ein wichtiges Detail in der Transformation, das oft vergessen wird. Bei der Untersuchung von Funktionswerten

f (t)

$f (t)$ Über das Intervall von 0 bis t wird die Funktion als periodisch betrachtet, wodurch Sie die Funktion zeitlich schmerzlos nach links oder rechts verschieben können. Daher haben wir das Recht, am Ende keinen neuen Wert einzufügen und den alten Wert von Anfang an zu löschen, sondern die Daten im Puffer zyklisch zu ersetzen.

Aus Gründen der Übersichtlichkeit können Sie in Tabellenform schreiben, wie sich der Puffer ändert:

t = 0	f (0)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 1	f (10)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 2	f (10)	f (11)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 3	f (10)	f (11)	f (12)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t = 4	f (10)	f (11)	f (12)	f (13)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)

Sie können schreiben, wie sich die zeitliche Transformation von t ₁ nach t ₂ ändert:

F_{t} = F_{t - 1} + D e l t a F D e l t a F : D e l t a c_{k} = f r a c 1 T i n t G r e n z e n_{t_{1}}^{t_{2}} (f_{t} (x) - f_{t - 1} (x)) e^{- 2 p i i k x / T} d x

$F_t = F_ {t-1} + \ Delta F \\ \ Delta F: \ Delta c_k = \ frac {1} {T} \ int \ Grenzen_ {t_1} ^ {t_2} (f_t (x) - f_ { t-1} (x)) e ^ {- 2 \ pi ikx / T} dx$

Wert

F_{t - 1} (x)

$F_ {t-1} (x)$ ist das Ergebnis der vorherigen Konvertierung und die Komplexität der Berechnung

D e l t a f (x)

$\ Delta f (x)$ hängt nicht von der Größe des Zeitfensters ab und ist daher konstant. Infolgedessen wird die Komplexität der Konvertierung sein

O (n)

$O (n)$ ^* weil Alles, was uns bleibt, ist, einmal durch das Frequenzfenster zu gehen und die Änderungen für die T-Samples anzuwenden, die sich im Laufe der Zeit geändert haben. Ich möchte auch Ihre Aufmerksamkeit auf die Tatsache lenken, dass die Chancen

e^{- 2 p i i k x / T}

$e ^ {- 2 \ pi ikx / T}$ kann im Voraus berechnet werden, was einen zusätzlichen Produktivitätsgewinn ergibt, und es verbleiben nur zwei Operationen im Zyklus: Subtrahieren von reellen Zahlen und Multiplizieren einer reellen Zahl mit einer komplexen, in der Praxis sind beide Operationen einfach und billig.

Um das Bild zu vervollständigen, muss nur noch der Ausgangszustand angegeben werden, aber hier ist alles einfach:

F_{0} (x) = 0

$F_0 (x) = 0$

* - Natürlich bleibt die endgültige Komplexität der gesamten Transformation so

O (n^{2})

$O (n ^ 2)$ Es wird jedoch schrittweise über n Iterationen ausgeführt, während der Puffer aktualisiert wird.

O (n)

$O (n)$ - Dies ist die Komplexität der Aktualisierung der Daten, aber genau das ist es, was wir brauchen (bei Verwendung der FFT die Komplexität jeder Transformation

O (n l o g (n))

$O (nlog (n))$ )

Aber was ist, wenn Sie tiefer graben. Oder die zweite n loswerden

Ich möchte sofort reservieren, dass die nächsten Schritte nur dann anwendbar sind, wenn Sie nicht vorhaben, die inverse Transformation für das Ergebnis durchzuführen (um das Signal zu korrigieren oder eine Impulsantwort zu erhalten). Zunächst möchte ich Sie daran erinnern, dass wir durch die Konvertierung ein symmetrisches Datenarray erhalten, mit dem wir die Anzahl der Konvertierungen sofort um die Hälfte reduzieren können.

Lassen Sie uns nun den resultierenden Datensatz unter den Bedingungen des Problems analysieren. Wir haben eine Reihe komplexer Zahlen, von denen jede die Amplitude und Phase der Schwingungen bei einer bestimmten Frequenz beschreibt. Die Häufigkeit kann durch die Formel bestimmt werden:

f [j] = j f r a c f d N

$f [j] = j \ frac {fd} {N}$ für

j < f r a c N 2

$j <\ frac {N} {2}$ . Lassen Sie uns den Schritt des Frequenzfensters für unsere Daten bewerten:

D e l t a f = f r a c f d N

$\ Delta f = \ frac {fd} {N}$ Für N = 2 ¹⁴ : 2,93 Hz (und für 2 ¹⁶ : 0,73 Hz). Somit erhalten wir im Bereich von 1 kHz bis 2 kHz 341 Ergebnisse. Versuchen Sie unabhängig zu bewerten, wie viele Daten für N = 65536 im Bereich von 8 kHz bis 16 kHz liegen. Sehr viel, richtig? Sehr viel! Benötigen wir so viele Daten? Bei den Problemen der Anzeige der Frequenzeigenschaften von Soundsystemen lautet die Antwort natürlich Nein. Andererseits ist für die Analyse im Niederfrequenzbereich ein kleiner Schritt sehr nützlich. Vergessen Sie nicht, dass noch ein Zeitplan für diese Bände vor uns liegt (

f r a c N 2

$\ frac {N} {2}$ ) in eine für Menschen lesbare Form konvertieren (Mittelwertbildung, Spline oder Glättung) und auf dem Bildschirm anzeigen. Und bei hohen Frequenzen, selbst wenn ein 4K-Bildschirm vorhanden ist und das Diagramm im Vollbildmodus mit der logarithmischen Frequenzachse angezeigt wird, beträgt die Schrittgröße schnell weniger als 1 Pixel.

Erfahrungsgemäß können Sie feststellen, dass es ausreicht, nur 48 Punkte pro Oktave zu haben und die Daten etwas flüssiger und gemitteler zu machen. Ich empfehle, bei 96 anzuhalten. Im Audiofrequenzbereich von 20 Hz bis 20 kHz ist es einfach, nur 10 Oktaven zu zählen: 20, 40, 80 , 160, 320, 640, 1280, 2560, 5120, 10240, 20480, von denen jedes in eine bestimmte Anzahl von Teilbändern unterteilt werden kann (vergessen Sie nicht, dass die Partition geometrisch und nicht arithmetisch erfolgen sollte), daher ist es mehr als ausreichend, die Konvertierung nur für durchzuführen 960 Frequenzen zu bekommen Performan dass in 16 ... 65-mal kleiner als die ursprüngliche Version.

Wenn wir beide Ansätze kombinieren, erhalten wir die konstante Komplexität des Datenaktualisierungsalgorithmus

O (1)

$O (1)$ .

Quadratischer Honig und ein Löffel Teer

Jetzt können wir das aus der Komplexität sicher sagen

O (n^{2})

$O (n ^ 2)$ Wir kamen zur Komplexität

O (1)

$O (1)$ mit zwei einfachen Life-Hacks:

Nach der Analyse des Problems stellten wir fest, dass die Daten schrittweise hinzugefügt werden und die Zeitspanne der vollständigen Aktualisierung des Zeitfensters viel höher ist als die Zeitspanne der Transformationen, und berechneten anschließend die Differenz der Fourier-Transformation.
vom arithmetischen Schritt im Frequenzfenster auf nur durch die angegebenen Werte begrenzt, was die Anzahl der Konvertierungen drastisch reduzieren kann.

Aber natürlich wäre das Leben wirklich ein Märchen, wenn nicht eines, aber. Die Anwendung dieser beiden Ansätze ermöglichte es uns, die CPU wirklich zu entladen, um zu erraten, dass sie die Fourier-Transformation berechnet und die Ergebnisse auch mit auf dem Bildschirm anzeigt

N = 2^{20}

$N = 2 ^ {20}$ es war schwierig. Die Bestrafung ließ sich jedoch nicht warten, wenn Ihre Signale in der Realität nicht periodisch sind (und dies ist erforderlich, um die richtigen Konvertierungsergebnisse zu erhalten) und es nicht möglich ist, die entsprechende Fenstergröße auszuwählen. Es wird erforderlich, verschiedene Fensterfunktionen zu verwenden, sodass Sie den ersten Schritt nicht mehr vollständig ausführen können. Die Praxis hat gezeigt, dass die Verwendung von Fensterfunktionen bei der Untersuchung von Signalen mit einer Frequenz von weniger als kritisch ist

0, 1 f_{d}

$0,1f_d$ . Bei hohen Frequenzen verringert die Anzahl der Perioden, die in das Zeitfenster fallen, die Verzerrungen, die sich aus dem Vorhandensein einer Lücke erster Ordnung (zwischen f (0) und f (N-1)) in der ursprünglichen Funktion ergeben, erheblich.

Am Ende lehnte ich auch den zweiten Schritt ab und kehrte zurück zur FFT, weil Der Gewinn bei dieser Aufgabe war bereits gering.

Abschließend

Der erste Ansatz kann angewendet werden, wenn Ihre Daten ausgesprochen periodisch sind und im Laufe der Zeit unter Verwendung eines großen Zeitfensters analysiert werden müssen, das, wie ich mich erinnere, nicht Grad 2 sein muss, d. H. jede natürliche Zahl.
Der zweite Ansatz ist anwendbar (auch unter Berücksichtigung von Fensterfunktionen), wenn nur ein bestimmter kleiner Satz von Frequenzen in den Daten analysiert wird.

Leider blieb dies für mich in diesem Problem nur eine kleine mathematische Unterhaltung, aber ich hoffe, dass es Sie dazu inspiriert, an Feiertagen andere Algorithmen im Hinblick auf Änderungen der Eingabedaten im Laufe der Zeit zu studieren :)

Literatur

Bild aus dem Manga von Michio Shibuya. “AUFREGENDE MATHEMATIK. Fourier-Analyse

Fourier-Transformation. Das schnelle und das wütende

N loswerden

Aber was ist, wenn Sie tiefer graben. Oder die zweite n loswerden

Quadratischer Honig und ein Löffel Teer

Abschließend

Literatur

More articles: