Calcul de Convolution en Numpy/Scipy
Profilage certains travaux informatiques, je suis en train de faire m'a montré qu'un goulot d'étranglement dans mon programme a une fonction qui fondamentalement fait (npnumpy,spscipy):
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)
les deux signaux ont une forme (C, N) où C est le nombre d'ensembles de données (habituellement moins de 20) et N est le nombre d'échantillons dans chaque ensemble (environ 5000). Le calcul pour chaque ensemble (ligne) est complètement indépendant de tout autre ensemble.
je me suis rendu compte que ceci c'était juste une simple convolution, j'ai donc essayé de le remplacer par:
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
...pour voir si j'ai eu les mêmes résultats. Mais je ne l'ai pas fait, et mes questions sont:
- Pourquoi pas?
- Est-il une meilleure manière de calculer l'équivalent de <!--10?
(je me rends compte que mix2 n'aurait probablement pas été plus rapide tel quel, mais cela aurait pu être un bon point de départ pour la parallélisation.)
voici le script complet que j'ai utilisé pour vérifier rapidement ceci:
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)
3 réponses
donc j'ai testé ceci et je peux maintenant confirmer quelques choses:
1) numpy.convolve n'est pas circulaire, ce qui est ce que le code fft vous donne:
2) le FFT n'est pas muni d'un tampon interne d'une puissance de 2. Comparer très différentes vitesses des opérations suivantes:
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 normalisation n'est pas une différence -- si vous faites une convolution circulaire naïve en additionnant a(k)*b(i-k), vous obtiendrez le résultat du code FFT.
La chose est rembourrage à une puissance de 2 va changer la réponse. J'ai entendu des histoires qu'il y a des façons de gérer cela en utilisant habilement les facteurs premiers de la longueur (mentionné mais pas codé dans les recettes numériques) mais je n'ai jamais vu les gens font réellement cela.
scipy.signal.fftconvolve convolve par FFT, c'est du code python. Vous pouvez étudier le code source et corriger la fonction mix1.
comme mentionné plus haut, le scipy.signal.la fonction convolve n'effectue pas de convolution circulaire. Si vous voulez une convolution circulaire exécutée dans l'espace réel (par opposition à l'utilisation de fft) je suggère d'utiliser le scipy.ndimage.de convolution de la fonction. Il a un paramètre de mode qui peut être réglé sur 'wrap' ce qui en fait une convolution circulaire.
for idx, row in enumerate(outputs):
outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode='wrap')