Die Faltung ist eine nützliche mathematische Operation, die in verschiedenen Gebieten für die Signal- und Bildverarbeitung, Statistik und mehr angewendet wird. Akustikingenieure verwenden die Faltung häufig zur Verarbeitung akustischer Signale, um die gewünschten Informationen zu erhalten oder den Klang besser zu interpretieren. In diesem Blog-Beitrag stellen wir drei Methoden zur Durchführung von Faltungen in der Software COMSOL Multiphysics® vor. Wir werden die Implementierung der Faltung anhand der Tiefpassfilterung einer Raumimpulsantwort diskutieren und ein Beispiel für die Auralisierung der Raumakustik zeigen.
Faltung: Definition und Methoden
Physikalisch gesehen gibt die Faltung Aufschluss über den Grad der Überlappung zwischen zwei Signalen, die im Zeitbereich, im Frequenzbereich oder im räumlichen Bereich dargestellt werden. Für Signale im Zeitbereich ist sie mathematisch definiert als
Dabei sind t und \tau die Zeitvariable bzw. die Dummy-Variable, die für die Zeitintegration verwendet werden, und ast ist der Faltungsoperator.
Zu jedem t berechnet die Gleichung das Zeitintegral des Produkts einer Funktion, f(\tau), in ihrer ursprünglichen Form mit der anderen Funktion, g(t-\tau), die von t gespiegelt und verschoben wird. Die Operation ist kommutativ, d.h. das Ergebnis bleibt dasselbe, egal welche Funktion gespiegelt und verschoben wird. Das Integral sollte für alle Werte der Verschiebung (t) durchgeführt werden, um das Faltungsergebnis als Funktion von t zu erhalten.
Diese Integralform, hier Faltungsintegral genannt, eignet sich für die Verarbeitung glatter und stetiger Funktionen. Für diskrete Daten, wie digitalisierte Tonsignale, erfordert diese Methode ein sehr anspruchsvolles numerisches Integrationsschema, was oft zu einem hohen Rechenaufwand führt. Um das zu vermeiden, kann die folgende diskrete Faltung für die Verarbeitung diskreter Signale verwendet werden:
Dabei ist \Delta t das Abtastintervall.
Durch Annäherung des Integrals an eine Summierung diskreter Abtastwerte kann die diskrete Faltung schneller berechnet werden als das Faltungsintegral.
Wenn die Fourier-Transformation der beiden Signale existiert, kann die diskrete Faltung mit Hilfe der Fast-Fourier-Transformation (FFT) auf der Grundlage des Faltungssatzes effizienter ausgewertet werden:
Dabei sind \mathcal{F} und \mathcal{F}^{-1} die Operatoren der Fourier-Transformation bzw. der inversen Fourier-Transformation. Das Faltungstheorem nutzt die Tatsache, dass die Fourier-Transformierte der Faltung zweier Funktionen im Zeitbereich äquivalent zum Produkt der Fourier-Transformierten der Signale (im Frequenzbereich) ist. Die moderne Echtzeit-Faltungstechnologie verwendet im Allgemeinen die Fast-Fourier-Transformation, weil sie sehr effizient ist.
Tiefpassfilterung der Raumimpulsantwort
In diesem Abschnitt wird gezeigt, wie die Faltung mit Hilfe des Faltungsintegrals, der diskreten Faltung und des Faltungssatzes in COMSOL Multiphysics® implementiert werden kann. Wir werden die wesentlichen Schritte dieser Methoden anhand eines Beispiels aus dem Tutorial-Modell Akustik eines kleinen Konzertsaals unter Anwendung der Tiefpassfilterung auf eine Impulsantwort (impulse response, IR) erläutern, da die Tiefpassfilterung häufig verwendet wird, um nur die tiefen Basstöne in Sounddesigns zu extrahieren. Diese IR-Daten werden geladen und in einer Interpolationsfunktion gespeichert (in diesem Fall Interpolation 1). Ein Plot der Daten ist unten abgebildet.
Ein Plot der Raumimpulsantwort.
Die Dauer und die Abtastfrequenz der Daten betragen 2 s bzw. 44.100 Hz.
Als Tiefpassfilter wird ein Butterworth-Filter 4. Ordnung verwendet. Die Übertragungsfunktion des Filters, TF(f), ist definiert als:
mit
wobei f und \omega die Frequenz bzw. die Winkelfrequenz darstellen. \omega_c ist die Winkelfrequenz bei der Grenzfrequenz. \Pi ist der Produktoperator, der das Produkt einer Sequenz darstellt.
Im Beispiel wird der Filter 4. Ordnung mit einer Grenzfrequenz von 2 kHz durch die Funktion Analytic definiert (in diesem Fall Analytic 1). Die folgenden Abbildungen zeigen, wie die Funktion definiert ist, den Frequenzgang der Funktion in Bezug auf ihren Real- und Imaginärteil und die Verstärkungskennlinien des durch eine solche Funktion definierten Filters.
Die Übertragungsfunktion des Butterworth-Filters 4. Ordnung wird mit Hilfe des Features Analytic definiert.
Der Frequenzgang der Übertragungsfunktion in Bezug auf ihren Real- und Imaginärteil.
Die Verstärkungskennlinie des Tiefpassfilters.
Die Verstärkungskennlinie weist darauf hin, dass der Tiefpassfilter bis zur Grenzfrequenz von 2 kHz einen flachen Durchlassbereich hat (an diesem Punkt dämpft der Filter die Eingangsleistung um die Hälfte oder 3 dB). Oberhalb der Grenzfrequenz nimmt die Filterverstärkung mit einer Rate von -24 dB pro Oktave ab.
Nun werden wir die drei Methoden zur Durchführung der Faltung in COMSOL Multiphysics® erläutern.
Methode 1: Direkte Behandlung des Faltungsintegrals
Wir beginnen mit der direkten Behandlung des Faltungsintegrals. Hier benötigt man zwei Signale aus dem Zeitbereich zur Eingabe. Das erste Signal, die Raumimpulsantwort, ist bereits als Interpolation 1 definiert worden. Das zweite Signal, der Tiefpassfilter, ist im Frequenzbereich definiert. Wir müssen es in den Zeitbereich umwandeln, indem wir eine inverse diskrete Fourier-Transformation (IDFT) des Signals durchführen. Das geschieht mit den folgenden Schritten.
Schritt 1
Erstellen Sie einen Grid 1D Datensatz und nennen Sie ihn Grid 1D_f. Das erzeugt ein Frequenzgitter, das in einem Frequenzbereich (-22050 Hz, 22050 Hz) angegeben ist, der dem Frequenzspektrum der Raumimpulsantwort-Daten entspricht. Es wird für den Plot der Signale im Frequenzbereich verwendet.
Einstellung für den Datensatz Grid 1D_f.
Schritt 2
Führen Sie die IDFT des Tiefpassfilters mit einem Plot Function im Feature 1D Plot Group mit dem Datensatz Grid 1D_f durch. Im Abschnitt Output des Einstellungsfensters wählen Sie Discrete Fourier transform von der Liste Display, dann wählen Sie Real part von der Liste Show und markieren Inverse transform.
Die Einstellungen und der Plot der IDFT des Tiefpassfilters.
Schritt 3
Fügen Sie Plot-Daten zu Table 1 hinzu.
Schritt 4
Definieren Sie eine Interpolationsfunktion (Interpolation 2) mit den in Tabelle 1 gespeicherten Daten.
Die Einstellung der Interpolation 2.
Jetzt haben wir zwei Signale im Zeitbereich, die für die Faltung verwendbar sind. Obwohl in der Definition von einem unendlichen Zeitintervall -\infty, \infty ausgegangen wird, muss das Integral nur in einem endlichen Zeitintervall von 0-2 s berechnet werden, da beide Eingangssignale außerhalb dieses Zeitbereichs auf Null gesetzt werden. Das Verfahren für das Faltungsintegral ist im Folgenden dargestellt.
Schritt 1
Definieren Sie im Knoten Geometry ein Intervall, dessen Anfangs- und Endwerte dem Integralintervall (0-2 s) entsprechen.
Erzeugen des Intervalls.
Schritt 2
Definieren Sie den Integrationsoperator intop1 auf dem Intervall.
Definition des Integraloperators.
Schritt 3
Diskretisieren Sie das Intervall mit einem einheitlichen Liniennetz, dessen Größe dem Abtastintervall der ursprünglichen Raumimpulsantwort-Daten entspricht.
Diskretisierung des Intervalls.
Schritt 4
Führen Sie eine Studie des Typs Stationary durch, damit die Interpolationsfunktionen und der Operator intop1 im Abschnitt Results verwendet werden können.
Schritt 5
Erstellen Sie einen Grid 1D Datensatz und nennen Sie ihn Grid 1D_t. Damit erzeugen Sie ein Zeitraster zur Definition des Zeitsignals im Bereich von 0-2 s mit einer Abtastfrequenz von 44.100 Hz.
Einstellungen für den Datensatz Grid 1D_t.
Schritt 6
Führen Sie das Faltungsintegral in einem 1D-Plot mit Grid 1D_t als Quelldatensatz aus.
Ausdruck für das Faltungsintegral unter Verwendung des Operators dest.
Dabei sind IR und Filter_IDFT die in Interpolation 1 und Interpolation 2 definierten Interpolationsfunktionen der Raumimpulsantwort und der Impulsantwort des Tiefpassfilters.
Es ist zu beachten, dass der Operator dest die Auswertung der Funktion am Zielpunkt statt am Quellpunkt erzwingt.
Methode 2: Standardmäßige diskrete Faltung
Die diskrete Faltung wird in der Praxis häufig verwendet. Um in COMSOL Multiphysics® eine diskrete Faltung durchzuführen, ist die Verwendung des Application Builders unerlässlich. Mit dem Method Editor im Application Builder können Sie Code von Java zum Erstellen von Methoden verwenden, die Sie ausführen können, um Operationen im Modellbaum zu automatisieren oder zu erweitern. In diesem Blog-Beitrag wird die Methode zur Durchführung einer Faltung mit in Tabellen gespeicherten Daten erläutert.
Das folgende Bild zeigt den implementierten Code und die Einstellungen und die folgende Tabelle enthält die im Code definierten Variablen.
Java-Methode zur Durchführung der diskreten Faltung mit in Tabellen gespeicherten Daten.
| Name | Typ | Beschreibung |
|---|---|---|
| a | Double (Array 2D) | Daten aus der ersten Eingabetabelle |
| b | Double (Array 2D) | Daten aus der zweiten Eingabetabelle |
| c | Double (Array 2D) | Ergebnis der Faltung |
| ndata1 | Integer (Skalar) | Länge von a |
| ndata2 | Integer (Skalar) | Länge von b |
| ndata | Integer (Skalar) | Länge von c |
| dt | Double (Skalar) | Abtastintervall |
| js | Integer (Skalar) | Startwert des Summationsindexes |
| je | Integer (Skalar) | Endwert des Summationsindexes |
Im Programm Java definierte Variablen.
Der Code führt die diskrete Faltung mit den in zwei Eingabetabellen gespeicherten Daten durch und exportiert das Ergebnis in eine Ausgabetabelle. Die folgenden Hinweise können zum besseren Verständnis des Codes verwendet werden:
- 2.-4. Zeilen: Die Daten und die Länge der Daten werden aus der ersten Eingabetabelle extrahiert.
- 6.-8. Zeilen: Die Daten und die Länge der Daten werden aus der zweiten Eingabetabelle extrahiert.
- 11.-13. Zeilen: Die Länge des Ergebnisses wird festgelegt, das Array mit doppelter Genauigkeit für die Speicherung des Faltungsergebnisses erstellt und das Abtastintervall definiert.
- 18.-28. Zeilen: Die For-Schleife wird ausgeführt, um die diskrete Faltung vom ersten bis zum letzten Zeitschritt zu starten.
- 30.-33. Zeilen: Das Ergebnis wird in die Ausgabetabelle mit der Bezeichnung Discrete convolution result exportiert.
Beachten Sie, dass die Länge der Ergebnisdaten die Summe der Länge der beiden Eingabetabellen abzüglich 1 ist und dass der Startwert (js) und der Endwert (je) des Summenindex so definiert sind, dass js kleiner als die Länge der zweiten Tabelle und je kleiner als die Länge der ersten Tabelle ist (siehe 22. bzw. 23. Zeile des Codes).
Um das Programm ausführen zu können, muss der Plot von zwei Zeitsignalen in Tabellen gespeichert werden. Die Plot-Daten der IDFT des Tiefpassfilters werden in Table 1 gespeichert. Für die Raumimpulsantwort muss im Feature 1D Plot Group mit dem Datensatz Grid 1D_t ein Funktionsplot der Raumimpulsantwort erstellt und die Plotdaten in Table 2 kopiert werden.
Die Tags der Tabellen können im Fenster Settings eines dem Modellbaum hinzugefügten Features Method Call eingegeben werden. Nachdem alle Einstellungen vorgenommen wurden, können Sie die diskrete Faltung durch Klicken auf die Schaltfläche Run im Einstellungsfenster von Method Call durchführen.
Eingabe von Tabellen-Tags im Fenster Settings des Features Method Call.
Methode 3: Diskrete Faltung auf der Grundlage des Faltungstheorems
Die letzte Methode führt die Faltung auf der Grundlage des Faltungstheorems, also unter Verwendung der Diskrete Fourier-Transformation, durch. In diesem Fall werden die DFT der Raumimpulsantwort und die Übertragungsfunktion des Tiefpassfilters verwendet. Die Methode ist wie folgt:
Schritt 1
Erstellen Sie sowohl für den Real- als auch für den Imaginärteil der DFT der Raumimpulsantwort jeweils einen Plot der Art Function in der gleichen 1D-Plot-Gruppe. Wählen Sie in den Einstellungen Discrete Fourier transform aus der Liste Display. Anschließend wählen Sie aus der Liste Show die Option Real part bzw. Imaginary part, um den Real- bzw. Imaginärteil der DFT der Raumimpulsantwort darzustellen. Die Auswahlfelder sind Standardeinstellungen.
DFT der Raumimpulsantwort in einem 1D-Plot.
Schritt 2
Kopieren Sie die Plot-Daten der Real- und Imaginärteile der DFT der Raumimpulsantwort in Table 3 bzw. Table 4.
Schritt 3
Definieren Sie die Interpolation 3 und die Interpolation 4 mit den in Table 3 bzw. Table 4 gespeicherten Daten.
Schritt 4
Führen Sie die Faltung mittels Berechnung der IDFT des punktweisen Produkts aus der DFT der Raumimpulsantwort und dem Tiefpassfilter durch. Die punktweise Multiplikation und die IDFT werden beide in einem einzigen Plot der Art Function unter Verwendung des Datensatzes Grid 1D_f durchgeführt.
Die IDFT des Produkts aus der DFT der Raumimpulsantwort und des Tiefpassfilters. Das Auswahlfeld Inverse transform im Fenster Settings ist markiert, um die inverse Transformation durchzuführen.
Dabei sind real_IR und imag_IR der Realteil und der Imaginärteil der DFT der Raumimpulsantwort und werden in Interpolation 3 bzw. Interpolation 4 definiert. Butterworth ist die in Analytic 1 definierte Übertragungsfunktion des Tiefpassfilters.
Beachten Sie, dass die Methode des Faltungstheorems eine zirkuläre Faltung durchführt, sodass eine kreisförmige Überlappung entsteht, wenn die Länge des Gitterintervalls nicht ausreicht.
Ergebnis der Tiefpassfilterung
Die folgenden Diagramme zeigen die mit den drei Methoden berechneten tiefpassgefilterten Raumimpulsantworten, aus denen eine gute Übereinstimmung zwischen den drei Methoden hervorgeht.
Ein Vergleich der Faltungsergebnisse, die mit dem Faltungsintegral, der diskreten Faltung und dem Faltungstheorem berechnet wurden.
Die Faltungsergebnisse liegen im Bereich von 0,02-0,07 s.
Die nächste Abbildung zeigt einen Vergleich der Spektren der Raumimpulsantwort vor und nach der Filterung. Es ist zu erkennen, dass das gefilterte Spektrum oberhalb von 2 kHz, der Grenzfrequenz des Tiefpassfilters, abnimmt. Dieses Ergebnis spricht für den Erfolg der Faltungsimplementierung.
Spektren der Raumimpulsantwort vor und nach der Filterung.
Anwendung auf die Auralisation
Kommen wir nun zum Beispiel der Auralisation. In diesem Beispiel geht es um die Faltung der Raumimpulsantwort und einer in einer Absorberkammer aufgezeichneten Tonaufnahme. Unter der Annahme eines linearen, zeitinvarianten Systems kann die Reaktion auf eine beliebige Eingabe durch Faltung von IR und Eingangssignal bewertet werden. Basierend auf dieser Theorie bewerten Akustiker Schallfelder auditiv durch die Faltung von mittels Rechenverfahren und reflexionsfreiem Schall simulierten IR. Dieser Simulationsprozess, welcher sich von der Erzeugung der Schalldaten bis zum Hören erstreckt, wird Auralisation genannt. Das Beispielmodell (erreichbar über die Schaltfläche im nächsten Abschnitt) führt die Auralisation eines Trompetenklangs im Modell des kleinen Konzertsaals mit der oben beschriebenen Standardmethode der diskreten Faltung durch. Sie können das Faltungsergebnis auch als WAV-Audiodatei zum Abspielen in einem Audio- oder Media-Player exportieren. Im Folgenden können Sie den reflexionsfreien und den hallenden Trompetenklang vergleichen.
Die Trompete, wie sie klingt, wenn sie in einer Absorberkammer gespielt wird.
Die Trompete, wie sie im kleinen Konzertsaal klingt, mit Auralisation.
Das obige Beispiel ist ein monauraler Klang, der sich von einem binauralen Klang unterscheidet, den wir normalerweise hören. Wir können jedoch realistischere Klänge erzeugen, indem wir die kopfbezogene Übertragungsfunktion oder die Techniken zur Schallfeldwiedergabe wie Ambisonics kombinieren (siehe Ref. 1). Solche Beispiele werden wir in einem zukünftigen Blogbeitrag vorstellen.
Nächste Schritte
Um das in diesem Blog-Beitrag besprochene Modell näher zu erkunden, klicken Sie auf die Schaltfläche unten, die Sie zur Application Gallery bringt.
Weitere Lektüre
In den folgenden Blog-Beiträgen erfahren Sie mehr über die Datenverarbeitung bei der Akustikmodellierung:
- Computing Transient Sound Pressure Levels in COMSOL®
- Modeling Room Acoustics Using a Hybrid Approach
Referenzen
- M. Vorländer, Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality, Springer Science & Business Media: Berlin/Heidelberg, 2007.
Der reflexionsfreie Schall wird bereitgestellt durch The Open AIR Library unter CC BY 4.0.
Oracle und Java sind eingetragene Marken von Oracle und/oder seinen Tochtergesellschaften.
Kommentare (0)