CODEWITHSUNDEEP
pythonnumpyscipysignal-processing
user2857354
user2857354

Reputation: 41

scipy.signal.fftconvolve doesn't give the required results

I have a question regarding python's fftconvolve. In my current research I've been required to calculate some convolution between two functions. To do so I'm calculating it using fourier transform (which I used numpy.fft and normalize it) . The thing is that if I want to compare it using fftconvolve package, it fails to give the correct results. Here is my code:

#!/usr/bin/python
import numpy as np
from scipy.signal import fftconvolve , convolve 

def FFT(array , sign):
  if sign==1:
    return np.fft.fftshift(np.fft.fft(np.fft.fftshift(array))) * dw / (2.0 * np.pi)
  elif sign==-1:
    return np.fft.fftshift(np.fft.ifft(np.fft.fftshift(array))) * dt * len(array)


def convolve_arrays(array1,array2,sign):
  sign = int(sign)
  temp1 = FFT(array1 , sign,)
  temp2 = FFT(array2 , sign,)
  temp3 = np.multiply(temp1 , temp2)
  return  FFT(temp3 , -1 * sign) / (2. * np.pi) 

""" EXAMPLE """ 

dt    = .1
N     = 2**17
t_max = N * dt / 2
time  = dt * np.arange(-N / 2 , N / 2 , 1)

dw    = 2. * np.pi / (N * dt)
w_max = N * dw / 2.
w     = dw * np.arange(-N / 2 , N / 2 , 1)

eta_fourier = 1e-10




Gamma   = 1.
epsilon = .5
omega   = .5


G    = zeros(N , complex)
G[:] = 1. / (w[:] - epsilon + 1j * eta_fourier)

D    = zeros(N , complex)
D[:] = 1. / (w[:] - omega + 1j * eta_fourier) - 1. / (w[:] + omega + 1j * eta_fourier)

H    = convolve_arrays(D , G , 1)     
J    = fftconvolve(D , G , mode = 'same') * np.pi  / (2. * N)

If you plot the real/imaginary part of H, J you'll see a shift in the w axes and also I had to multiply the J's results in order to get somehow close (but still not) to the correct results.

Any suggestions?

Thanks!

Upvotes: 4

Views: 4175

Answers (1)

Andrew
Andrew

Reputation: 2892

Boundary conditions are important when you compute convolutions.

When you convolve two signals, the edges of the result depend on what values you assume outside the edges of the inputs. fftconvolve computes convolution assuming zero-padded boundaries.

Take a look at the source code of fftconvolve. Notice the shenanigans they go through to achieve zero-padded boundary conditions, in particular, these lines:

size = s1 + s2 - 1

...

fsize = 2 ** np.ceil(np.log2(size)).astype(int) #For speed; often suboptimal!
fslice = tuple([slice(0, int(sz)) for sz in size])

...

ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()

...

return _centered(ret, s1) #strips off padding

This is good stuff! It's probably worth reading fftconvolve's code carefully, and a good education if you want to understand Fourier-based convolution.

Brief sketch

The forward FFT zero-pads each signal to prevent periodic boundary conditions:

a = np.array([3, 4, 5])
b = np.fft.ifftn(np.fft.fftn(a, (5,))).real
print b #[ 3.  4.  5.  0.  0.]

the inverse FFT of the product of the forward FFTs gives a padded result:

a = np.array([3, 4, 5])
b = np.array([0., 0.9, 0.1])
b = np.fft.ifftn(np.fft.fftn(a, (5,))*
                 np.fft.fftn(b, (5,))
                 ).real
print b #[ 0.   2.7  3.9  4.9  0.5]

and the _centered function strips off the extra padding pixels at the end (assuming you use the mode='same' option).

Upvotes: 7

Related Questions