CODEWITHSUNDEEP
pythonnumpysignal-processingfftconvolution
Jack Rolph
Jack Rolph

Reputation: 595

NumPy Convolution Theorem

I am new to convolution and would therefore like to prove convolution theorem to myself by convolving two 1D signals together using FFT. However, my code is not consistent with the result obtained from np.convolve.

I have used the naive implementation from this solution:

import numpy as np

def Convolution(array,kernel):
    return np.real(np.fft.ifft( np.fft.fft(array)*np.fft.fft(kernel) ))

 a_flat = [ 1., 2., 3., 0., 0., 4., 5., 6., 0., 0. , 7.,  8.,  9.]
 k_flat = [ 1,2 ,1 ,0,0 ,0 ,0,0 ,0 ,0,-1,-2,-1]

 my_convolution = Convolution(a_flat, k_flat)
 np_convolution = np.convolve(a_flat, k_flat)

 print(my_convolution)
 print("") 
 print(np_convolution)

The output reads:

[ 19.  10.   4.  -5. -17. -13.   7.  13.  -5. -26. -20.   9.  24.]

[  1.   4.   8.   8.   3.   4.  13.  20.  17.   6.   6.  18.  24.  18.   6.  -4. -13. -20. -17.  -6.  -7. -22. -32. -26.  -9.]

I am clearly missing something. Can someone point me to my oversight?

Upvotes: 1

Views: 1036

Answers (3)

hotpaw2
hotpaw2

Reputation: 70743

The result of a convolution is N+M-1, e.g. longer than either of the inputs. So your FFT needs to be that long or longer.

An FFT/IFFT will wrap the fast convolution result around, and mix it up into a circular convolution. But if you pad the data with lots of zeros on the end(s), the mix will be easy to unmix.

Upvotes: 1

roadrunner66
roadrunner66

Reputation: 7941

What @hotpaw2 said. Always better to graph it: 

import numpy as np
import matplotlib.pyplot as p
%matplotlib inline

def Convolution(array,kernel):
    return np.real(np.fft.ifft( np.fft.fft(array)*np.fft.fft(kernel) ))

a_flat = [ 1., 2., 3., 0., 0., 4., 5., 6., 0., 0. , 7.,  8.,  9.]
k_flat = [ 1,2,1,0,0,0,0,0,0,0,-1,-2,-1]

a_flat= np.pad(a_flat, (25, 25), 'constant', constant_values=(0, 0)).tolist()
k_flat= np.pad(k_flat, (25, 25), 'constant', constant_values=(0, 0)).tolist()


my_convolution = Convolution(a_flat, k_flat)
np_convolution = np.convolve(a_flat, k_flat)

fig,ax = p.subplots(3,figsize=(12,5))     

ax[0].plot(a_flat)
ax[1].plot(k_flat)
ax[2].plot(np.roll(my_convolution, 30),'.-',lw=0.5,label='myconv');  # arbitrary shift here
ax[2].plot(np.roll(np_convolution,  0),'.-',lw=3,alpha=0.3,label='npconv');
p.legend()

Nice example, btw.

enter image description here

Upvotes: 2

code-lukas
code-lukas

Reputation: 1651

Be aware of where your anchor is in your kernel, usually it's the median, in your case the 0 in the middle, you should still ensure it's correct. Referring to @Paul R 's answer, make up your mind what kind of padding you want to use (zero padding, replicate the border regions, etc.) as this will influence your output.

Upvotes: 1

Related Questions