# convolution with gaussian kernel using fft

46 views (last 30 days)
LM on 13 Jun 2020
Edited: David Goodmanson on 15 Jun 2020
Hey,
I'm really no pro in Matlab so I've got a few difficulties with the following task.
I'm trying to implement diffusion of a circle through convolution with the 2d gaussian kernel. The convolution is between the Gaussian kernel an the function u, which helps describe the circle by being +1 inside the circle and -1 outside. The Gaussian kernel is . I've tried not to use fftshift but to do the shift by hand. Also I know that the Fourier transform of the Gaussian is with coefficients depending on the length of the interval.
But with my code, there happens no diffusion at all and I don't understand what I've done wrong. I would relly appreciate some help.
I have the follwoing code with meshsize nxn:
function gaussian(n)
length = 1; %length of the interval
x = (length/n)*(0:n-1);
[X1,X2] = meshgrid(x,x); %grid
K = [0:n/2-1,-n/2:-1]; [K1,K2] = meshgrid(K,K); %fftshift by hand
A = K1.^2 + K2.^2; %coefficients for the Fourier transform of the Gaussian kernel
dt = 0.01;
R0 = 0.4; %radius of the circle
%initial condition
u = u_0(X1,X2,n,length,R0);
show(X1,X2,u,length);
for i=0:dt:1
Gaussian = (length/n)^2*exp(-dt*A); %pre-factor of the discrete fourier transform
convolution = sign(real((length/n)^(-2)*ifft2( (length/n)^2*fft2(u) .* Gaussian ))); %here I'm solving the convolution with fft2
show(X1,X2,convolution,length); drawnow
end
function show(X1,X2,u,length)
surf(X1,X2,u);
shading flat; axis([0 length 0 length -1 1 -1 1]);
%inital condition
function val = u_0(X1,X2,n,length,r)
u = zeros(n,n);
A = (X1-length/2).^2 + (X2-length/2).^2;
for i=1:n
for j=1:n
if A(i,j) < r^2
u(i,j) = 1;
else
u(i,j) = -1;
end
end
end
val = u;

Image Analyst on 13 Jun 2020
Why aren't you using conv2()?
LM on 13 Jun 2020
It is an application to the theory of fast fourier transform, that's why I'm using fft.

David Goodmanson on 14 Jun 2020
Hi LM
The code below takes your approach but modifies some of the details.  The array sizes are odd x odd since you get the greatest amount of symmetry that way.  As you can see, the code uses a lot of fftshift and ifftshift. For odd arrays, fftshift is not its own inverse and neither is ifftshift its own inverse. But they are inverses of each other. All you need to remember is that fftshift takes k=0 or x=0 to the center of the array, and ifftshift takes them both down to the origin so that fft and ifft work properly.  the k array needs a factor of 2pi compared to what you have, since k = 2pi (1/lambda), analogous to w = 2pi f.
 I used a function that is 0 outside the circle instead of -1. Using -1 gives a huge peak at k =0 in the k domain, but using 0 gives a much better plot of what is going on in the k domain. You can change to -1 at the indicated point in the code and the result is good for that case.  You don't have to worry about normalization of the kernel in k space. That's because the function is exp(-k.^2*t). You know the function should = 1 when t = 0, because there should be no effect. So the constant in front has to be 1.
n = 81; % should be odd
xymax = 10; % sets the spacial domain
t = 2;
nov2 = (n-1)/2;
xx = (-nov2:nov2)*(xymax/nov2);
kk = (-nov2:nov2)*(nov2/xymax)*2*pi/n;
[x y] = meshgrid(xx,xx);
[k q] = meshgrid(kk,kk);
u0 = 0*ones(size(x)); % change to -1 for original problem
u0((x.^2+y.^2)<=1) = 1; % 1 inside the disc
G = exp(-(k.^2+q.^2)*t); % gaussian kernel in k space
u0f = fftshift(fft2(ifftshift(u0))); % fourier transform of u0
ucheck = fftshift(ifft2(ifftshift(u0f))); % transform back as a check
max(max(abs(u0 -ucheck))) % should be small
ut = fftshift(ifft2(ifftshift(u0f.*G))); % evolution in time
figure(1)
surf(x,y,u0); view([0 0 1])
figure(2)
surf(k,q,abs(u0f))
figure(3)
surf(x,y,ucheck); view([0 0 1])
figure(4)
surf(x,y,ut); view([0 0 1])

LM on 14 Jun 2020
Thanks a lot, now it's working! I just don't really understand your : why should k needs a factor 2pi and what du you mean by w?
David Goodmanson on 14 Jun 2020
In the time domain, (analogous to the space domain), if the initial array is time then the fft produces a function of frequency. And you can calculate a frequency array for the output. But w (i.e. omega, the circular frequency) = 2*pi*f, So to get an array in w instead of f, you multiply the frequencies by 2*pi. The analogous situations for time and space are:
time t f w = 2*pi*f
space x (1/lambda) k = 2*pi*(1/lambda)
so obtaining k requires a factor of 2*pi.
LM on 15 Jun 2020