0
$\begingroup$

I use the pretty simple example used in many books to understand the convolution in the frequency domain. But I cannot find the real results.

import numpy as np

img = np.array([[0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0]])
kernel = np.array([[0, 0, 0, 0, 0],
                   [0, 1, 2, 3, 0],
                   [0, 4, 5, 6, 0],
                   [0, 7, 8, 9, 0],
                   [0, 0, 0, 0, 0]])

fft_img = np.fft.fft2(img)
fft_kernel = np.fft.fft2(kernel)

fshift = (fft_img* fft_kernel)  # Convelution theory at frequency domain

img_back = np.fft.ifft2(fshift)

img_back = np.real(img_back)
img_back = img_back.astype(np.int)

print(img_back)

The output of the img_back is

[[9 0 0 7 8]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [3 0 0 0 1]
 [6 0 0 4 5]]

However, the answer should be

[[0 0 0 0 0]
 [0 9 8 7 0]
 [0 6 5 4 0]
 [0 3 2 1 0]
 [0 0 0 0 0]]

So, is there someone who can help me? I don't know the reason why. And how to come true the ideal result?

$\endgroup$

1 Answer 1

1
$\begingroup$

This is because the identity for the convolution is

[[1, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0]]

You get the correct result by replacing

img_back = np.real(img_back)

by

img_back = np.fft.ifftshift(np.round(np.real(img_back)))

It is important to np.round (to the nearest integer) before casting to np.int, otherwise you will may have wrong results due to rounding errors e.g. np.array([0.999999999999]).astype(np.int) = 0

Notes about fftshift

The the fftshift and ifftshift operations are equivalent to a convolution, we could express fftshift(X) as $S * X$. Convolutions are commutative and associative $a * (b * c) = (a * b) * c$, and $a * b = b * a$, it is immediate in the frequency from the multiplication operator properties.

Lets say $X$ is your image $H$ is your kernel

you can check that $S^{-1} * X$, or $X = S$, is precisely the term that I mentioned in the beginning of the answer is the identity for the convolution

What you have from my solution above is $S^{-1} * (H * X)$, the solution you metioned in the comments is $H * (S^{-1} * X)$, and there would be one other variation $(S^{-1} * H) * X$ that would represents calling fftshift on the kernel instead of the image. I already saw some doing ffthisft in both inputs and then ifftshift to the output.

Alternative forms...

$(S^{-1} * X) * H = X * (S^{-1} * H) = S * (S^{-1} * X) * (S^{-1} * H)$

$\endgroup$
8
  • $\begingroup$ Thank you very much. I found when I insert the code img = np.fft.ifftshift(img) before fft_kernel = np.fft.fft2(img), this also can get the same result. But I still don't understand what you said about the impulse identify for the convolution. $\endgroup$
    – Guo Ray
    Feb 26, 2021 at 8:02
  • $\begingroup$ updated the answer explaining this case $\endgroup$
    – Bob
    Feb 26, 2021 at 13:41
  • $\begingroup$ Thank you very much. It is really clear. But I still don't understand what you said about the impulse identify for the convolution, or I should ask why we need an ifftshift(X) for this operation. You know there is not fftshift(X). $\endgroup$
    – Guo Ray
    Mar 2, 2021 at 1:03
  • $\begingroup$ If thinking about circular shifting of negative indices is not helping, think about two signals starting at with duration N/2, centered at N/2, it means they have non-zero values from N/4 to 3N/4. once you convolve them the result will be possibly non-zero in the range N/2 to 3N/2, but you compute the FFT using only N samples, you assign the interval N/2 to 3N/2, to the indices 0 to N. This is why we use fftshift. $\endgroup$
    – Bob
    Mar 2, 2021 at 7:40
  • $\begingroup$ Thake a look at this answer in the 1D context, it includes plots that may help to build your intuition dsp.stackexchange.com/questions/73471/… $\endgroup$
    – Bob
    Mar 2, 2021 at 7:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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