# Implementing Convolution in Frequency Domain?

Suppose, we have a bitmap image represented as a 2D integer array, int [,] image2D; whose FFT is Complex[,] fftImage2D;

Suppose, we have an kernel represented as a 2D integer array, int [,] kernel2D; whose FFT is Complex[,] fftKernel2D;

We know that, the convolution (in spatial domain) of image2D and kernel2D would be,

int Rows = image2D.GetLength(0);
int Cols = image2D.GetLength(1);

for(int i=0 ; i<Rows ; i++)
{
for(int j=0 ; j<Cols ; j++)
{
//sweep the kernel2D across image2D
//...........................
}
}


We also know that, convolution in frequency domain would be, multiplication between fftImage2D and fftKernel2D.

How can I do this multiplication?

How can I multiply two Complex [,] type 2D arrays of different dimensions? I have understood the theory. My problem is practical implementation. As I described in the question,

1. Are DFT of the image and DFT of the kernel going to be of different sizes? I guess so. So, how can I multiply them element by element?

2. In my code, each of the DFTs are represented by 2D Complex numbers. Should, I multiply them according to complex-number's multiplication rule? Probably yes. But, only when their dimensions are same. Right?

IDFT the smaller or both of the DFTs if needed. Zero pad one or both of the kernel and image to make them the same dimension and size. Re-DFT as needed, and now you can complex multiply the 2 DFT arrays element-by-element because they will now be the same size.

• But, if I zero-pad the smaller one (in most cases, the kernel), should I keep the kernel at the center and surround them with zeros, or should I shift the kernel in one of the corners? – user18425 Jun 12 '16 at 22:23
• Depends on how you want the result offset. I typically zero-pad both up to a factorable-into-small-primes size that is equal or greater than N+M-1, with the originals centered in the zero padding, or fftshift one or both. – hotpaw2 Jun 12 '16 at 23:34
• What is the common/standard/general procedure? – user18425 Jun 12 '16 at 23:38

This (kind of) question has been possibly asked and answered a few times, but for your convenience let me put another.

Given a 2D discrete-space sequence (representing a sampled image) of $$x[n_1,n_2]$$ of length $$N_1$$ along $$n_1$$-axis from $$n_1=0$$ to $$n_1=N_1-1$$ and length $$N_2$$ along $$n_2$$-axis from $$n_2=0$$ to $$n_2=N_2-1$$, and a 2D discrete-space filter's impulse response sequence $$h[n_1,n_2]$$ of length $$L_1$$ along $$n_1$$-axis and length $$L_2$$ along $$n_2$$-axis similarly.

Then, their 2D discrete-space convolution is defined as

$$y[n_1 , n_2] = x[n_1 , n_2] * h[n_1 , n_2]$$ $$y[n_1 , n_2] = \sum_{m_1} \sum_{m_2} {x[m_1,m_2] h[n_1-m_1 , n_2-m_2]}$$

and has the lengths $$M_1 = N_1 + L_1 -1$$ along $$n_1$$ axis starting from $$n_1 = 0$$ and $$M_2 = N_2 + L_2 - 1$$ along $$n_2$$ axis similary.

By remembering the well known convolution - multiplication theorem of the DSFT (discrete space Fourier transform) one can also compute the result $$y[n_1,n_2]$$ in the frequency domain as :

$$x[n_1,n_2] \star h[n_1,n_2] \longleftrightarrow X(e^{j\omega_1}e^{j\omega_2}) H(e^{j\omega_2}e^{j\omega_2})$$

The theorem used above computes a linear convolution between the input sequences through multiplying their theoretical DSFTs. However its computer implementation is based on the practical DFT sequences $$X[k_1,k_2]$$ and $$H[k_1,k_2]$$. Since their multiplication performs an equivalent circular convolution of $$x[n_1,n_2]$$ and $$h[n_1,n_2]$$, one should choose the DFT lengths carefully to avoid aliasing in the result :

$$K_1 \geq M_1 = N_1 + L_1 -1$$ $$K_2 \geq M_2 = N_2 + L_2 -1$$

Eventhough any $$K_1$$, $$K_2$$ satisfying above will do, it's logical to set them either to $$M_1$$ and $$M_2$$ or the next higher power of $$2$$ for performance reasons. The three step process begins by computing the two $$M_1 \times M_2$$ - point 2D-DFTs, $$X[k_1,k_2]$$ and $$H[k_1,k_2]$$ of the input sequences $$x[n_1,n_2]$$ and $$h[n_1,n_2]$$. Then the resulting $$M_1 \times M_2$$ point 2D DFTs are multiplied sample by sample to produce the the output 2D-DFT sequence $$Y[k_1,k_2]$$. Finally the multiplicand $$Y[k_1,k_2]$$ is transformed back into space domain by an $$M_1 \times M_2$$ point inverse DFT to produce the expected resulting sequence $$y[n_1,n_2]$$ in a much shorter processing time when the filter's impulse response sequence length is longer than some decidable number.

Implementing the above in matlab/octave can be done with the following line:

Iy = real( ifft2( fft2(Ix, M2,M1).*fft2(h,M2,M1) , M2, M1) );

• I think you need to shift the data. Otherwise you're doing a Convolution of shifted kernel. Remember the [0, 0] of the Filter is actually in the middle yet MATLAB treat as it was in the first sample ([1, 1] in this case). – Royi Apr 28 '17 at 11:38
• @Royi Have no idea of what you mean by the shift of data, but the above single MATLAB line is all that's required to perform the frequency domain equivalent of the time domain convolution between two signals with appropriate sizes. There is no explicit index referrence in the algorithm. Both signals are assumed to begin from the index [0,0] in their mathematical description and therefore from index [1,1] of the their respective MATLAB array counterparts. The fact that MATLAB indexing is beginning from n=1 does not create a problem here, unlike in cases where you explicitly use the indices. – Fat32 Apr 28 '17 at 13:57
• When some builds a 2D Filter, usually, the axis center is in the middle of the Filter. For instance, imagine Gaussian Kernel, the center of the filter is in its center element of the matrix. Yet, you algorithm above in MATLAB will yield Convolution with respect to the top left element. Hence your result will be effectively a shifted version of the desired result. – Royi Apr 28 '17 at 14:36