Note that there are some explanatory texts on larger screens.

plurals
  1. POscipy.signal.fftconvolve doesn't give the required results
    text
    copied!<p>I have a question regarding python's <code>fftconvolve</code>. 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 <code>numpy.fft</code> and normalize it) . The thing is that if I want to compare it using <code>fftconvolve</code> package, it fails to give the correct results. Here is my code:</p> <pre><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) </code></pre> <p>If you plot the real/imaginary part of <code>H</code>, <code>J</code> you'll see a shift in the <code>w</code> axes and also I had to multiply the <code>J</code>'s results in order to get somehow close (but still not) to the correct results.</p> <p>Any suggestions?</p> <p>Thanks! </p>
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload