Second derivative using fft

Question

All, I am trying to take the laplacian of the following function:

g(x,y) = 1/2cx^2+1/2dy2

The laplacian is c + d, which is a constant. Using FFT I should get the same ( in my FFT example I am padding the function to avoid edge effects).

Here is my code:

Define a 2D function
n = 30 # number of points
Lx = 30 # extension in x
Ly = 30 # extension in x
dx = n/Lx # Step in x 
dy = n/Ly # Step in x 
c=4
d=4
x=np.arange(-Lx/2,Lx/2)
y=np.arange(-Ly/2,Ly/2)
g = np.zeros((Lx,Ly))
lapg = np.zeros((Lx,Ly))

for j in range(Ly):
    for i in range(Lx):
        g[i,j] = (1/2)*c*x[i]**2 + (1/2)*d*y[j]**2
        lapg[i,j] = c + d

kxpad = 2*np.pi*np.fft.fftfreq(2*Lx,d=dx)
#kxpad = (2*np.pi/(2*Lx))*np.arange(-2*Lx/2,2*Lx/2)
#kxpad = np.fft.fftshift(kxpad)

#kypad = (2*np.pi/(2*Ly))*np.arange(-2*Ly/2,2*Ly/2)
#kypad = np.fft.fftshift(kypad)

kypad = 2*np.pi*np.fft.fftfreq(2*Ly,d=dy)
kpad = np.zeros((2*Lx,2*Ly))


for j in range(2*Ly):
    for i in range(2*Lx):
        kpad[i,j] =  math.sqrt(kxpad[i]**2+kypad[j]**2)

kpad = np.fft.fftshift(kpad)

gpad = np.zeros((2*Lx,2*Ly))


gpad[:Lx,:Ly] = g # Filling main part of g in gpad

gpad[:Lx,Ly:] = g[:,-1::-1] # Filling the last 3 columns of gpad with g flipped

gpad[Lx:,:Ly] = g[-1::-1,:]# Filling the last 3 lines of gpad with g flipped

gpad[Lx:,Ly:] = g[-1::-1, -1::-1]# Filling the last 3 lines and last 3 columns of gpad with g flipped in line and column

rdFFT2D = np.zeros((Lx,Ly))

gpadhat = np.fft.fft2(gpad)
dgpadhat = -(kpad**2)*gpadhat #taking the derivative iwFFT(f)
rdpadFFT2D = np.real(np.fft.ifft2(dgpadhat))
rdFFT2D = rdpadFFT2D[:Lx,:Ly]

[

First image is the plot of the original function g(x,y), 2nd image is the analytical laplacian of g and 3rd image is the sugar loaf in Rio de Janeiro( lol ), actually it is the laplacian using FFT. What Am I doing wrong here?

Edit : Commenting on ripple effect. Cris you mean the ripple effect due to the set_zlimit in the image below?Just to remember you that the result should be 8.

Edit 2 : Using non-symmetrical x and y values, produce the two images.

Answer 1

The padding will not change the boundary condition: You are padding by replicating the function, mirrored, four times. The function is symmetric, so the mirroring doesn't change it. Thus, your padding simply repeats the function four times. The convolution through the DFT (which you're attempting to implement) uses a periodic boundary condition, and thus already sees the input function as periodic. Replicating the function will not improve the convolution results at the edges.

To improve the result at the edges, you would need to implement a different boundary condition, the most effective one (since the input is analytical anyway) is to simply extend your domain and then crop it after applying the convolution. This introduces a boundary extension where the image is padded by seeing more data outside the original domain. It is an ideal boundary extension suitable for an ideal case where we don't have to deal with real-world data.

This implements the Laplace though the DFT with greatly simplified code, where we ignore any boundary extension, as well as the sample spacing (basically setting dx=1 and dy=1):

import numpy as np
import matplotlib.pyplot as pp

n = 30 # number of points
c = 4
d = 4
x = np.arange(-n//2,n//2)
y = np.arange(-n//2,n//2)
g = (1/2)*c*x[None,:]**2 + (1/2)*d*y[:,None]**2

kx = 2 * np.pi * np.fft.fftfreq(n)
ky = 2 * np.pi * np.fft.fftfreq(n)
lapg = np.real(np.fft.ifft2(np.fft.fft2(g) * (-kx[None, :]**2 - ky[:, None]**2)))

fig = pp.figure()
ax = fig.add_subplot(121, projection='3d')
ax.plot_surface(x[None,:], y[:,None], g)
ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(x[None,:], y[:,None], lapg)
pp.show()

---

Edit: Boundary extension would work as follows:

import numpy as np
import matplotlib.pyplot as pp

n_true = 30 # number of pixels we want to compute
n_boundary = 15 # number of pixels to extend the image in all directions
c = 4
d = 4

# First compute g and lapg including boundary extenstion
n = n_true + n_boundary * 2
x = np.arange(-n//2,n//2)
y = np.arange(-n//2,n//2)
g = (1/2)*c*x[None,:]**2 + (1/2)*d*y[:,None]**2
kx = 2 * np.pi * np.fft.fftfreq(n)
ky = 2 * np.pi * np.fft.fftfreq(n)
lapg = np.real(np.fft.ifft2(np.fft.fft2(g) * (-kx[None, :]**2 - ky[:, None]**2)))

# Now crop the two images to our desired size
x = x[n_boundary:-n_boundary]
y = y[n_boundary:-n_boundary]
g = g[n_boundary:-n_boundary, n_boundary:-n_boundary]
lapg = lapg[n_boundary:-n_boundary, n_boundary:-n_boundary]

# Display
fig = pp.figure()
ax = fig.add_subplot(121, projection='3d')
ax.plot_surface(x[None,:], y[:,None], g)
ax.set_zlim(0, 800)
ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(x[None,:], y[:,None], lapg)
ax.set_zlim(0, 800)
pp.show()

Note that I'm scaling the z-axes of the two plots in the same way to not enhance the effects of the boundary too much. Fourier-domain filtering like this is typically much more sensitive to edge effects than spatial-domain (or temporal-domain) filtering because the filter has an infinitely-long impulse response. If you leave out the set_zlim command, you'll see a ripple effect in the otherwise flat lapg image. The ripples are very small, but no matter how small, they'll look huge on a completely flat function because they'll stretch from the bottom to the top of the plot. The equal set_zlim in the two plots just puts this noise in proportion.