Breve sintesi: Come si calcola rapidamente la convoluzione finita di due array?Manufatti di Riemann sum in scipy.signal.convolve
campioni discreti Descrizione del problema
Sto provando avere la convoluzione finita di due funzioni f (x), g (x) definito da
Per ottenere questo, ho preso delle funzioni e trasformati in campi di lunghezza steps:
xarray = [x * i/steps for i in range(steps)]
farray = [f(x) for x in xarray]
garray = [g(x) for x in xarray]
ho quindi cercato di calcolare il conv oluzione usando la funzione scipy.signal.convolve. Questa funzione fornisce gli stessi risultati dell'algoritmo conv suggerito here. Tuttavia, i risultati differiscono notevolmente dalle soluzioni analitiche. La modifica dell'algoritmo conv per utilizzare la regola trapezoidale fornisce i risultati desiderati.
Per illustrare questo, ho lasciato
f(x) = exp(-x)
g(x) = 2 * exp(-2 * x)
i risultati sono:
Qui Riemann rappresenta una semplice somma di Riemann, trapezoidal è una versione modificata dell'algoritmo di Riemann di utilizzare la regola trapezoidale, scipy.signal.convolve è la funzione scipy e analytical è la convoluzione analitica.
Ora diamo g(x) = x^2 * exp(-x) ed i risultati diventano:
Qui 'ratio' è il rapporto tra i valori ottenuti SciPy ai valori analitici. Quanto sopra dimostra che il problema non può essere risolto rinormalizzando l'integrale.
La domanda
E 'possibile utilizzare la velocità di SciPy ma mantenere i migliori risultati di una regola trapezoidale o devo scrivere un estensione C per raggiungere i risultati desiderati?
Un esempio
Basta copiare e incollare il codice sottostante per vedere il problema che sto incontrando. I due risultati possono essere portati ad un accordo più stretto aumentando la variabile steps. Credo che il problema sia dovuto a artefatti tratti dalle somme di Riemann di destra, perché l'integrale è sovrastimato quando aumenta e si avvicina di nuovo alla soluzione analitica mentre diminuisce.
EDIT: Ora ho inserito l'algoritmo originale 2 come un confronto che dà gli stessi risultati come la funzione scipy.signal.convolve.
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()
Grazie per il vostro tempo!
che potrebbe essere utile se si fornisce un [ completa l'esempio minimo] (http: // sscce.org /) che riproduce il problema, per escludere errori banali come una divisione intera viene utilizzato dove deve essere utilizzata la vera divisione. – jfs
se si sposta l'array 'sci' a destra (di un passo) e si normalizza quindi le soluzioni [hanno un aspetto simile] (http://i403.photobucket.com/albums/pp111/uber_ulrich/convolve.png): [' sci = np.r_ [0, sci [: steps-1]] * 0.86'] (https://gist.github.com/ac45fb5d117d8ecb66a3). Sembra che ci siano diverse definizioni di cosa 'convolve()' significa (fft, circolare), vedere la discussione in [Convolution computations in Numpy/Scipy] (http://stackoverflow.com/q/6855169/4279). – jfs
Grazie per i riferimenti al thread convoluzione. Il turno giusto sembra un modo interessante per risolvere il problema. Puoi dirmi come hai ottenuto il fattore di normalizzazione di 0,86? Ho incluso l'algoritmo di somma di Riemann originale per illustrare il motivo per cui ritengo che questo sia un artefatto numerico piuttosto che una diversa definizione di cosa significhi convoluzione. –