How to multiply the spectra of two images of different sizes?

Question

This is not a programming issue. But I am sure that this is something that is widely known and understood in this community.

I have an image, x and a much smaller image, y, and I need to convolve them by multiplying their FFT. But since they are not the same size, I do not know how to perform multiplication in the frequency domain.

I take the (two-dimensional) FFT x (which is an integer matrix of dimensions 4096 x 4096), which gives me a representation of the frequency domain X (which is a matrix of complex numbers, and I think its dimension is 2048 x 2048).

Similarly, I take (a two-dimensional FFT y (which is an integer matrix of dimension 64 x 64), which gives me a representation of the frequency domain Y (which is also a matrix of complex numbers, and I think this size is 32 x 32).

I use the fourn function in Numerical Recipes, so my input matrices x and y must be collapsed into one-dimensional arrays, which are replaced by their discrete Fourier transforms X and Y. The fact is that although this is a two-dimensional image problem, I work with one-dimensional arrays.

If I were trying to drill two images of the same size, x and y. Everything would be very simple:

X = FFT(x) Y = FFT(y) Z = X * Y (term by term multiplication) Convolution of x and y = IFFT(Z) 

But if X and Y are different lengths, how do I do the multiplication?

One possibility is to skip y in order to have the same dimensions as x. But that seems terribly ineffective. Another possibility is to align Y to have the same dimensions as X. But I don't know what that means in frequency space.

Here's another way to ask this question: if I want to drill two images of very different measurements using FFT so that I can multiply their spectra (represent the frequency domain), how do I multiply this?

Thanks,

~ Michael.

Answer 1

Filling a smaller array (convolution kernel, y in your case) with zeros according to the size of the input image (your matrix x) is the standard approach. It would be terribly inefficient if you were doing convolution in the spatial domain, but if you multiply the FFT, this is necessary, and the cost of calculating the FFT of the filled array is not so bad.

Answer 2

You are right in thinking that the intervals between the two frequencies should be the same. Take 1D example (I use Matlab syntax):

 N = 4096; M = 64; x = randn(N, 1); y = hann(M, 'symmetric'); zLinear = conv(x,y); zCircular = ifft( fft(x) .* fft(y,N) ); disp(max(abs(zLinear(65:4096) - zCircular(65:4096)))); 

The difference between the two methods is ~ 2e-14, so the rounding error. Note that you need to skip the first 64 patterns due to the difference between linear and circular convolution.

When calculating zCircular, pay attention to fft (y, N), which is the Matlab syntax for filling the signal y with zeros to N before calculating fft. This may be considered memory inefficient, but compare the speed:

linear convolution: 4096 times / add 64 each = 262144 times / add

circular convolution: 2 FFTs from 4096 + 1 complex multiplications by 2 * 4096 elements + 1 inverse FFT
= 3 * 4096 * log2 (4096) + 4096 * 6 = 172032 (provided that there are 6 operations for complex multiplication)

Basically, the NFTN FFT speed, even if you need three of them, is superior to N * M if M is very short.

EDIT Add Speed ​​Score for 2D Case

It’s worth adding that for 2D data, the speed advantage increases. 2D FFT performs N * N * log2 (N * N) operations; therefore, multiplying 3 FFTs + complex N ^ 2 by N = 4096 is 1.3 e10 operations. But direct convolution is N ^ 2 * M ^ 2 = 6.9e10 operations, about 50 times slower.