CODEWITHSUNDEEP
matlabsignal-processingfftoctavegaussian
koshy george
koshy george

Reputation: 680

Weird behavior when performing 2D convolution by the FFT

I am trying to construct the product of the FFT of a 2D box function and the FFT of a 2D Gaussian function. After, I find the inverse FFT and I am expecting a convolution of the two functions as a result. However, I am getting a weird one-sided result as seen below. The result appears on the bottom right of the subplot.

weird result

The Octave code I wrote to reproduce the above subplot as well as the calculations I performed to construct the convolution is shown below. Can anyone tell me what I'm doing wrong?

clear all;
clc;
close all;
% domain on each side is 0-9
L = 10;
% num subdivisions
N = 32;
delta=L/N;
sigma = 0.5;
% get the domain ready
[x,y] = meshgrid((0:N-1)*delta);
% since domain ranges from 0-(N-1) on both sdes
% we need to take the average
xAvg = sum(x(1, :))/length(x(1,:));
yAvg =  sum(y(:, 1))/length(x(:,1));
% gaussian
gssn =  exp(- ((x - xAvg) .^ 2 + (y - yAvg) .^ 2) ./ (2*sigma^2));


function ret = boxImpulse(a,b)
    n = 32;
    L=10;
    delta = L/n;
    nL = ((n-1)/2-3)*delta;
    nU = ((n-1)/2+3)*delta;
    if ((a >= nL) && (a <= nU) && ( b >= nL) && (b <= nU) )
        ret=1;
    else
        ret=0;
    end
    ret;
endfunction

boxResponse = arrayfun(@boxImpulse, x, y);

subplot(2,2,1);mesh(x,y,gssn); title("gaussian fun");
subplot(2,2,2);mesh(x,y, abs(fft2(gssn)) .^2); title("fft of gaussian");
subplot(2,2,3);mesh(x,y,boxResponse); title("box fun");
inv_of_product_of_ffts =  abs(ifft2(fft2(boxResponse) * fft2(gssn))) .^2 ;
subplot(2,2,4);mesh(x,y,inv_of_product_of_ffts); title("inv of product of fft");

Upvotes: 2

Views: 1243

Answers (1)

rayryeng
rayryeng

Reputation: 104464

Let's first address your most obvious errors. This is the way you are computing the convolution in frequency domain:

inv_of_product_of_ffts =  abs(ifft2(fft2(boxResponse) * fft2(gssn))) .^2 ;

The first problem is that you are using *, which is matrix multiplication. In frequency domain, element-wise multiplication is the equivalent to convolution in the frequency domain, so you need to use .* instead. The second problem is you also don't need the .^2 term and the abs operation. You're performing convolution, then finding the absolute value and squaring each term. This isn't required. Remove the abs and .^2 operations.

Therefore, what you need is:

inv_of_product_of_ffts = ifft2(fft2(boxResponse).*fft2(gssn)));

However, what you're going to get is this result. Let's place this in a new figure instead of a subplot*:

figure;
mesh(x,y,ifft2(fft2(boxResponse).*fft2(gssn))); 
title('Convolution... not right though');

enter image description here

You can see that it's the right result... but it's not centred.... why is that? This is actually one of the most common problems when computing convolution in the frequency domain. In fact, even the most experienced encounter this problem because they don't understand the internals of how the FFT works.

This is a consequence with how MATLAB and Octave compute the 2D FFT. Specifically, MATLAB and Octave defines a meshgrid of coordinates that go from 0,1,...M-1 for the horizontal and 0,1,...N-1 for the vertical, giving us a M x N result. This means that the origin / DC component is at the top-left corner of the matrix, not the centre as we usually define things. Specifically, the traditional 2D FFT defines the coordinates from -(M-1)/2, ..., (M-1)/2 for the horizontal and -(N-1)/2, ..., (N-1)/2 for the vertical.

Referencing the aforementioned, you defined your signal with the centre assuming that it's the origin and not the top-left corner. To compensate for this, you need to add an fftshift (MATLAB doc, Octave doc) so that the output of the FFT is now centred at the origin and not the top-left corner to bring things back to the way they were, and therefore you really need:

figure;
mesh(x,y,fftshift(ifft2(fft2(boxResponse).*fft2(gssn)))); 
title('Convolution... now it is right');

enter image description here

If you want to double check that we have the right result, you can perform direct convolution in the spatial domain and we can compare the results between the method in frequency domain.

This is how you'd compute the result directly in spatial domain:

figure;
mesh(x,y,conv2(gssn,boxResponse,'same'));
title('Convolution... spatial domain');

conv2 (MATLAB doc, Octave doc) performs 2D convolution between two signals, and the 'same' flag ensures that the output size is the largest of the two signals, which is either the Gaussian or Box filter. You'll see that it's the same curve, and I won't show it here for brevity.

However, we can compare both of the results and see if they're the same element-wise. A method to do this would be to determine if subtracting each element in the result is less than some threshold... say.. 1e-10:

>> out1 = conv2(boxResponse, gssn, 'same');
>> out2 = fftshift(ifft2(fft2(boxResponse).*fft2(gssn)));
>> all(abs(out1(:)-out2(:)) < 1e-10)

ans =

     1

This means that indeed both of these are the same.

Hope that makes sense! Remember, since you're defining your signal where the origin is at the centre, once you find the inverse FFT, you must shift things back so that the origin is now at the centre, not the top-left corner.

*: Minor note - All plots produced in this answer were generated with MATLAB R2015a. However, the code fully works in Octave - tested with Octave 4.0.0.

Upvotes: 7

Related Questions