Artefactos de la suma de Riemann en scipy.signal.convolve

Question

Breve resumen: ¿Cómo calculo rápidamente la convolución finita de dos matrices?

muestras discretas Descripción del problema

Estoy tratando de obtener la convolución finito de dos funciones f (x), g (x) definida por

Para lograr esto, que he tomado de las funciones y los convirtió en matrices de longitud steps:

xarray = [x * i/steps for i in range(steps)] 
farray = [f(x) for x in xarray] 
garray = [g(x) for x in xarray] 

luego trató de calcular el conv solución usando la función scipy.signal.convolve. Esta función proporciona los mismos resultados que el algoritmo conv sugerido here. Sin embargo, los resultados difieren considerablemente de las soluciones analíticas. La modificación del algoritmo conv para usar la regla trapezoidal da los resultados deseados.

Para ilustrar esto, vamos

f(x) = exp(-x) 
g(x) = 2 * exp(-2 * x) 

los resultados son los siguientes:

Aquí Riemann representa un simple suma de Riemann, trapezoidal es una versión modificada del algoritmo de Riemann a utilizar el regla trapezoidal, scipy.signal.convolve es la función scipy y analytical es la convolución analítica.

Ahora vamos g(x) = x^2 * exp(-x) y los resultados se convierten en:

Aquí 'relación' es la proporción de los valores obtenidos de scipy a los valores analíticos. Lo anterior demuestra que el problema no se puede resolver renormalizando la integral.

La pregunta

¿Es posible utilizar la velocidad de scipy pero conservan los mejores resultados de una regla trapezoidal o tengo que escribir una extensión C para lograr los resultados deseados?

Un ejemplo

Simplemente copia y pega el código siguiente para ver el problema que estoy encontrando. Los dos resultados pueden ponerse de acuerdo para aumentar la variable steps. Creo que el problema se debe a los artefactos de la mano derecha de las sumas de Riemann porque la integral se sobreestima cuando aumenta y se acerca a la solución analítica nuevamente a medida que disminuye.

EDITAR: Ahora he incluido el algoritmo original 2 como una comparación que da los mismos resultados que la función 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() 

¡Gracias por su tiempo!

Answer 1

Respuesta corta: Write it in C!

Respuesta larga

Usando el libro de cocina sobre numpy arrays Reescribí el método de convolución trapezoidal en C. Para utilizar el código C uno requiere tres archivos (https://gist.github.com/1626919)

• El código C (performancemodule. do).
• El archivo de instalación para compilar el código y hacerlo invocable desde python (performancemodulesetup.py).
• El archivo de Python que hace uso de la extensión C (performancetest.py)

El código debe ejecutarse a partir de la descarga de la siguiente manera

• Ajuste la ruta de inclusión en performancemodule.c.

• Ejecutar el siguiente

  pitón performancemodulesetup.py construir pitón performancetest.py

Es posible que tenga que copiar el archivo de biblioteca o performancemodule.soperformancemodule.dll en el mismo directorio que performancetest.py.

resultados y el rendimiento

Los resultados concuerdan perfectamente uno con el otro como se muestra a continuación:

El rendimiento del método C es incluso mejor que el método de convolución de scipy. Correr 10k circunvoluciones con longitud de la matriz 50 requiere

convolve (seconds, microseconds) 81 349969 
scipy.signal.convolve (seconds, microseconds) 1 962599 
convolve in C (seconds, microseconds) 0 87024 

Por lo tanto, la aplicación C es aproximadamente 1000 veces más rápido que la aplicación pitón y un poco más de 20 veces más rápido que la aplicación scipy (la verdad, la aplicación scipy es más versátil).

EDIT: Esto no resuelve la pregunta original exactamente pero es suficiente para mis propósitos.

Answer 2

o, para aquellos que prefieren numpy a C. Será más lento que la implementación C, pero son solo unas pocas líneas.

>>> t = np.linspace(0, maxtime-dt, 50) 
>>> fx = np.exp(-np.array(t)) 
>>> gx = 2*np.exp(-2*np.array(t)) 
>>> analytical = 2 * np.exp(-2 * t) * (-1 + np.exp(t)) 

esto parece trapezoidal en este caso (pero no el registro de la matemáticas)

>>> s2a = signal.convolve(fx[1:], gx, 'full')*dt 
>>> s2b = signal.convolve(fx, gx[1:], 'full')*dt 
>>> s = (s2a+s2b)/2 
>>> s[:10] 
array([ 0.17235682, 0.29706872, 0.38433313, 0.44235042, 0.47770012, 
     0.49564748, 0.50039326, 0.49527721, 0.48294359, 0.46547582]) 
>>> analytical[:10] 
array([ 0.  , 0.17221333, 0.29682141, 0.38401317, 0.44198216, 
     0.47730244, 0.49523485, 0.49997668, 0.49486489, 0.48254154]) 

más grande error absoluto:

>>> np.max(np.abs(s[:len(analytical)-1] - analytical[1:])) 
0.00041657780840698155 
>>> np.argmax(np.abs(s[:len(analytical)-1] - analytical[1:])) 
6