Artefakte von Riemann sum in scipy.signal.convolve
Kurze Zusammenfassung : Wie berechne ich schnell die endliche Faltung zweier Arrays?
Problembeschreibung
Ich versuche, die endliche Faltung von zwei Funktionen f (x), g (x) zu erhalten, die durch
definiert sind 
Um dies zu erreichen, habe ich diskrete Beispiele der Funktionen genommen und sie in Arrays der Länge steps umgewandelt:
Ich habe dann versucht, die Faltung mit der Funktion scipy.signal.convolve zu berechnen. Diese Funktion liefert die gleichen Ergebnisse wie der Algorithmus conv hier vorgeschlagen . Die Ergebnisse unterscheiden sich jedoch erheblich von analytischen Lösungen. Das Ändern des Algorithmus conv zur Verwendung der Trapezregel gibt die gewünschten Ergebnisse.
Um das zu verdeutlichen, lasse ich
%Vor%Die Ergebnisse sind:
Hier stellt Riemann eine einfache Riemann-Summe dar, trapezoidal ist eine modifizierte Version des Riemann-Algorithmus zur Verwendung der Trapezregel, scipy.signal.convolve ist die scipy-Funktion und analytical ist die analytische Faltung.
Lassen Sie nun g(x) = x^2 * exp(-x) und die Ergebnisse werden:
Hier ist "Verhältnis" das Verhältnis der Werte, die von scipy zu den analytischen Werten erhalten werden. Das Obige zeigt, dass das Problem nicht durch Renormierung des Integrals gelöst werden kann.
Die Frage
Ist es möglich, die Geschwindigkeit von scipy zu verwenden, aber die besseren Ergebnisse einer Trapezregel beizubehalten oder muss ich eine C-Erweiterung schreiben, um die gewünschten Ergebnisse zu erzielen?
Ein Beispiel
Kopieren Sie einfach den folgenden Code und fügen Sie ihn ein, um das Problem zu sehen, das mir begegnet. Die beiden Ergebnisse können durch Erhöhen der steps -Variable zu einer besseren Übereinstimmung gebracht werden. Ich glaube, dass das Problem auf Artefakte von rechts zurückzuführen ist, Riemann summiert sich, weil das Integral im Anstieg überschätzt wird und sich der analytischen Lösung wieder nähert, wenn sie abnimmt.
BEARBEITEN : Ich habe jetzt den ursprünglichen Algorithmus 2 als ein Vergleich, der die gleichen Ergebnisse wie die Funktion scipy.signal.convolve liefert.
Danke für Ihre Zeit!
2 Antworten
Kurze Antwort : Schreiben Sie es in C!
Lange Antwort
Mit dem Kochbuch über numpy arrays habe ich die Trapezfaltung in C umgeschrieben. Um die C Code eins benötigt drei Dateien ( Ссылка )
- Der C-Code (performancemodule.c).
- Die Setup-Datei, um den Code zu erstellen und von python (performancemodulsetup.py) aufrufbar zu machen.
- Die Python-Datei, die die C-Erweiterung (performancet.py) verwendet
Der Code sollte beim Herunterladen ausgeführt werden, indem Sie Folgendes tun
- Passen Sie den Include-Pfad in
performancemodule.c. an
-
Führen Sie Folgendes durch
python performancemodulsetup.py build python performancetest.py
Sie müssen möglicherweise die Bibliotheksdatei performancemodule.so oder performancemodule.dll in dasselbe Verzeichnis kopieren wie performancetest.py .
Ergebnisse und Leistung
Die Ergebnisse stimmen wie folgt überein:
Die Leistung der C-Methode ist sogar besser als die Convolve-Methode von Scipy. Das Ausführen von 10k-Faltungen mit Array-Länge 50 erfordert
%Vor%Somit ist die C-Implementierung ungefähr 1000 mal schneller als die python-Implementierung und ein bisschen mehr als 20 mal so schnell wie die scipy-Implementierung (zugegebenermaßen ist die scipy-Implementierung vielseitiger).
BEARBEITEN : Das löst die ursprüngliche Frage nicht genau, aber ist für meine Zwecke ausreichend.
oder, für diejenigen, die numpy C vorziehen. Es wird langsamer sein als die C-Implementierung, aber es sind nur ein paar Zeilen.
%Vor%sieht in diesem Fall wie trapezförmig aus (aber ich habe die Mathematik nicht überprüft)
%Vor%größter absoluter Fehler:
%Vor%Tags und Links python scipy convolution