2

I'm trying to perform linear convolutions in Python by comparing the results from FFTs and convolution functions.

Python's scipy.signal.fftconvolve automatically does the necessary zero padding. If we do the calculation using only FFTs, we add a length of zeros after our input signal. Is there a way to trim the output array of the FFT result? (The fftconvolve has a mode='same' option where the output array is the same size as the input.)

import numpy as np
from numpy.fft import fft,ifft
from scipy.signal import fftconvolve
import matplotlib.pyplot as plt

tf = 50
ti = -tf
N = 2**10
t = np.linspace(ti,tf,N)
T = (tf-ti)/N

x = np.exp(-4*t**2)*np.cos(t)
y = (t-2)/((t-2)**2+3**2)

z1 = ifft(fft(x,len(x)+len(y)-1)*fft(y,len(x)+len(y)-1)) # need zero-padding
z2 = fftconvolve(x,y,mode='same') # automatic zero-padding

plt.figure(0)
plt.plot(t,z2)
#plt.figure(1)
#plt.plot(t,z1)

Note: Another issue is if we change the outputs to:

z1 = ifft(fft(x,len(x)+len(y)-1)*fft(y,len(x)+len(y)-1))*T
z2 = fftconvolve(x,y)*T

(where T is a normalisation constant) then plot either z1 or z2, we get strange boundary effects:

enter image description here

This effect is explained in this DSP post, but I'm not sure how to apply the solutions to my system.

Community
  • 1
  • 1
Medulla Oblongata
  • 3,771
  • 8
  • 36
  • 75
  • For what reason do you do the zero padding? – Chiel Jul 01 '16 at 09:04
  • If I'm not mistaken, Python does circular convolution, so we need the zero padding to do a linear convolution. See also: http://dsp.stackexchange.com/questions/741/why-should-i-zero-pad-a-signal-before-taking-the-fourier-transform – Medulla Oblongata Jul 01 '16 at 09:07
  • @MedullaOblongata I can't give you a solution to your problem, but when I used convolution product during my Astrophysics internship, I used astropy and his convolve function : http://docs.astropy.org/en/stable/convolution/ – Essex Jul 01 '16 at 09:18

0 Answers0