Kısa özet: İki dizinin sonlu konvolüsyonunu hızlı bir şekilde nasıl hesaplayabilirim?Riemann'dan elde edilen eserler scipy.signal.convolve
iki fonksiyon f (x), g sonlu kıvrım elde etmek için çalışıyorum
(x) Bu amaca ulaşmak için
ile tanımlanan bir sorun açıklama ben almış farklı örnekler
xarray = [x * i/steps for i in range(steps)]
farray = [f(x) for x in xarray]
garray = [g(x) for x in xarray]
sonra dönüşümler hesaplamaya çalıştı: fonksiyonların ve uzunluğu steps dizilerine çevirmişlerdir scipy.signal.convolve işlevini kullanarak olution. Bu işlev, conv önerilen here algoritmasıyla aynı sonuçları verir. Bununla birlikte, sonuçlar analitik çözümlerden önemli ölçüde farklıdır. Trapezoidal kuralı kullanmak için conv algoritmasını değiştirmek istenen sonuçları verir.
Riemann basit Riemann toplamını temsil, trapezoidal kullanmak Riemann algoritmasının değiştirilmiş bir versiyonudur:
f(x) = exp(-x)
g(x) = 2 * exp(-2 * x)
izin trapezoidal kural, scipy.signal.convolve scipy işlevi ve analytical analitik kesişimidir.
Şimdi g(x) = x^2 * exp(-x) izin vermedi ve sonuç olmak:
'oranı' analitik değerlere SciPy elde edilen değerlerin oranıdır. Yukarıdakiler, sorunun integrali yeniden düzenleyerek çözülemeyeceğini göstermektedir.
soru
o scipy hızını kullanmak ama bir yamuk kuralı iyi sonuçlar korumak veya ben istenen sonuçları elde etmek C uzantı yazmak zorunda mümkün mü?
bir örnek
Sadece kopyalamak ve ben karşılaşmamdır sorunu görmek için aşağıdaki kodu yapıştırın. İki sonuç, steps değişkenini artırarak daha yakın bir anlaşmaya varılabilir. Problemin, sağ el Riemann toplamlarından esinlendiğine inanıyorum çünkü integral, gittikçe artarken, analitik çözümü tekrar azalırken tekrar yaklaşıyor.
DÜZENLEME: Şimdi scipy.signal.convolve fonksiyonu olarak aynı sonuçları veren bir karşılaştırma olarak orijinal algoritma 2 dahil ettik.
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import math
def convolveoriginal(x, y):
'''
The original algorithm from http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
'''
P, Q, N = len(x), len(y), len(x) + len(y) - 1
z = []
for k in range(N):
t, lower, upper = 0, max(0, k - (Q - 1)), min(P - 1, k)
for i in range(lower, upper + 1):
t = t + x[i] * y[k - i]
z.append(t)
return np.array(z) #Modified to include conversion to numpy array
def convolve(y1, y2, dx = None):
'''
Compute the finite convolution of two signals of equal length.
@param y1: First signal.
@param y2: Second signal.
@param dx: [optional] Integration step width.
@note: Based on the algorithm at http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
'''
P = len(y1) #Determine the length of the signal
z = [] #Create a list of convolution values
for k in range(P):
t = 0
lower = max(0, k - (P - 1))
upper = min(P - 1, k)
for i in range(lower, upper):
t += (y1[i] * y2[k - i] + y1[i + 1] * y2[k - (i + 1)])/2
z.append(t)
z = np.array(z) #Convert to a numpy array
if dx != None: #Is a step width specified?
z *= dx
return z
steps = 50 #Number of integration steps
maxtime = 5 #Maximum time
dt = float(maxtime)/steps #Obtain the width of a time step
time = [dt * i for i in range (steps)] #Create an array of times
exp1 = [math.exp(-t) for t in time] #Create an array of function values
exp2 = [2 * math.exp(-2 * t) for t in time]
#Calculate the analytical expression
analytical = [2 * math.exp(-2 * t) * (-1 + math.exp(t)) for t in time]
#Calculate the trapezoidal convolution
trapezoidal = convolve(exp1, exp2, dt)
#Calculate the scipy convolution
sci = signal.convolve(exp1, exp2, mode = 'full')
#Slice the first half to obtain the causal convolution and multiply by dt
#to account for the step width
sci = sci[0:steps] * dt
#Calculate the convolution using the original Riemann sum algorithm
riemann = convolveoriginal(exp1, exp2)
riemann = riemann[0:steps] * dt
#Plot
plt.plot(time, analytical, label = 'analytical')
plt.plot(time, trapezoidal, 'o', label = 'trapezoidal')
plt.plot(time, riemann, 'o', label = 'Riemann')
plt.plot(time, sci, '.', label = 'scipy.signal.convolve')
plt.legend()
plt.show()
Zaman ayırdığınız için teşekkür ederiz!
sağladığınız eğer yardımcı olabilecek bir [ en az örnek tamamla] (http: // sscce.org /) sorunu yeniden üreten, tam bölünme gibi önemsiz hataları dışlamak için doğru bölüm kullanılmalıdır. – jfs
Eğer sci dizisini sağa (tek adım) kaydırırsanız ve normalleştirirseniz o zaman çözümler [benzer görünebilir] (http://i403.photobucket.com/albums/pp111/uber_ulrich/convolve.png): [' sci = np.r_ [0, sci [: adım-1]] * 0.86 '] (https://gist.github.com/ac45fb5d117d8ecb66a3). “Convolve()” ın anlamı (fft, dairesel), [Numpy/Scipy'deki Convolution hesaplamaları] konusundaki tartışmaya (http://stackoverflow.com/q/6855169/4279) ilişkin farklı tanımlar vardır. – jfs
Evrişim iş parçacığı referansları için teşekkür ederiz. Doğru vites, problemi düzeltmenin ilginç bir yolu gibi görünüyor. 0.86 normalleşme faktörünü nasıl elde ettiğinizi söyleyebilir misiniz? Orijinal Riemann toplam algoritmasını, bunun, konvolüsyonun ne anlama geldiğinin farklı bir tanımından ziyade niçin sayısal bir artefakt olduğuna inandığımı göstermek için dahil ettim. –