3
$\begingroup$

I have the following convolution as part of a numerical simulation.

$$T(r)=\int \mathrm{d}^3r_2\, p(r_2)f(r_2)\alpha(r-r_2)\, .$$

My problem is that the analytical expressions for $f$ and $p$ do exist but, I have the expression for $\alpha$ only in the Fourier domain in the form of $\alpha(k)$. I planned to evaluate in the following way:

  1. Construct a grid using the mesh grid of $100\times100\times100$ using meshgrid and linspace in numpy
ran = linspace(-1,1,N_r)
x,y,z = meshgrid(ran,ran,ran) #position space
  1. Construct the components xf, yf, zf in the Fourier domain from x, y, z
xf = fftn(x)
yf = fftn(y)
zf = fftn(z)
  1. Find the Fourier transform of $f(r)\times p(r)$ using fftn in numpy
  2. Multiply it with $\alpha(k)$
  3. Taking the inverse Fourier transform using ifftn in numpy.

I am not very sure that the above method works and I actually failed to verify it properly. I tried using scipy.ndimage.convolve to compare the results with the inverse Fourier transform of the product in the Fourier domain. Is it correct as to what I am doing with the code? And is there a way where I can verify that a method is working with a simpler example?

Trying to verify:

I have tried the following to test the theory. Seems like it does not work. I expect that the result res_1 and res_2 to be the same. I also used the function real to truncate the tiny imaginary part which results from the fftn and ifftn functions.

x = linspace(-1,1,10)
xf = fftn(x)

def f(x):
    return x**2+x**3*sin(x)

def g(k):
    return k**2+k**3/(3-k**2)

g_k = g(xf)
g_x = real(ifftn(g_k))

res_1 = img_con(g_x,f(x))

res_2 = real(ifftn(g(xf)*fftn(f(x))))

print(res_1)
print(res_2)

Am I doing something which is wrong?

$\endgroup$
2
$\begingroup$

You need to pay attention that unless properly padded the Multiplication in the Frequency Domain (DFT) applies Circular Convolution while you're after Linear Convolution.

For practical examples and more information have a look on my answers:

  1. Kernel Convolution in Frequency Domain - Cyclic Padding.
  2. Question About Linear and Circular Convolution - 1D and 2D.
$\endgroup$
  • $\begingroup$ But, I do not want the size of the convolution to change. I want the convolution to give me a result with the same dimensions as f(x). That is not possible, here as this is a discrete convolution, and this problem does not arise if I am doing a continuous convolution right? $\endgroup$ – lattitude Aug 29 at 11:30
  • 1
    $\begingroup$ In discrete world, in order to achieve Linear Convolution in Time / Spatial Domain using multiplication in Frequency Domain one must apply padding to signals. There is no way around it. $\endgroup$ – Royi Aug 29 at 11:55
  • $\begingroup$ Then, how can I get the convolution to be of the same size as the matrix f(x) $\endgroup$ – lattitude Aug 30 at 4:26
  • $\begingroup$ Have a look on my answers on dsp.stackexchange.com/a/56025/128 and dsp.stackexchange.com/a/56031/128. $\endgroup$ – Royi Aug 30 at 6:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.