python - one - Cálculos de convolución en Numpy/Scipy

finite convolution

Question

Hacer un perfil del trabajo computacional que estoy haciendo me mostró que un cuello de botella en mi programa era una función que básicamente hizo esto ( np es numpy , sp es scipy ):

def mix1(signal1, signal2):
    spec1 = np.fft.fft(signal1, axis=1)
    spec2 = np.fft.fft(signal2, axis=1)
    return np.fft.ifft(spec1*spec2, axis=1)

Ambas señales tienen forma (C, N) donde C es el número de conjuntos de datos (generalmente menos de 20) y N es el número de muestras en cada conjunto (alrededor de 5000). El cálculo para cada conjunto (fila) es completamente independiente de cualquier otro conjunto.

Pensé que esto era solo una simple convolución, así que intenté reemplazarlo con:

def mix2(signal1, signal2):
    outputs = np.empty_like(signal1)

    for idx, row in enumerate(outputs):
        outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode=''same'')

    return outputs

... solo para ver si tengo los mismos resultados. Pero no lo hice, y mis preguntas son:

1. Por qué no?
2. ¿Hay una mejor manera de calcular el equivalente de mix1() ?

(Me doy cuenta de que mix2 probablemente no habría sido más rápido como está, pero podría haber sido un buen punto de partida para la paralelización).

Aquí está el script completo que usé para revisar rápidamente esto:

import numpy as np
import scipy as sp
import scipy.signal

N = 4680
C = 6

def mix1(signal1, signal2):
    spec1 = np.fft.fft(signal1, axis=1)
    spec2 = np.fft.fft(signal2, axis=1)
    return np.fft.ifft(spec1*spec2, axis=1)

def mix2(signal1, signal2):
    outputs = np.empty_like(signal1)

    for idx, row in enumerate(outputs):
        outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode=''same'')

    return outputs

def test(num, chans):
    sig1 = np.random.randn(chans, num)
    sig2 = np.random.randn(chans, num)
    res1 = mix1(sig1, sig2)
    res2 = mix2(sig1, sig2)

    np.testing.assert_almost_equal(res1, res2)

if __name__ == "__main__":
    np.random.seed(0x1234ABCD)
    test(N, C)

Answer 1

Así que probé esto y ahora puedo confirmar algunas cosas:

1) numpy.convolve no es circular, que es lo que te da el código fft:

2) FFT no se adapta internamente a una potencia de 2. Compare las velocidades muy diferentes de las siguientes operaciones:

x1 = np.random.uniform(size=2**17-1)
x2 = np.random.uniform(size=2**17)

np.fft.fft(x1)
np.fft.fft(x2)

3) La normalización no es una diferencia: si realiza una convolución circular ingenua sumando a (k) * b (ik), obtendrá el resultado del código FFT.

La cosa es que el relleno a una potencia de 2 va a cambiar la respuesta. He escuchado historias de que hay maneras de lidiar con esto utilizando inteligentemente los factores primos de la longitud (mencionados pero no codificados en Recetas Numéricas) pero nunca he visto a personas que realmente lo hagan.

Answer 2

Como se mencionó anteriormente, la función scipy.signal.convolve no realiza una convolución circular. Si desea que se realice una convolución circular en el espacio real (en contraste con el uso de fft), sugiero usar la función scipy.ndimage.convolve. Tiene un parámetro de modo que se puede configurar para "envolver" y convertirlo en una convolución circular.

for idx, row in enumerate(outputs):
    outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode=''wrap'')

Answer 3

scipy.signal.fftconvolve se convierte por FFT, es el código de Python. Puede estudiar el código fuente y corregir su función mix1.