46 views (last 30 days)

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;

David Goodmanson
on 14 Jun 2020

Hi LM

The code below takes your approach but modifies some of the details. [1] The array sizes are odd x odd since you get the greatest amount of symmetry that way. [2] 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. [3] the k array needs a factor of 2pi compared to what you have, since k = 2pi (1/lambda), analogous to w = 2pi f.

[4] 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. [5] 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])

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.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Direct link to this comment

https://ww2.mathworks.cn/matlabcentral/answers/547593-convolution-with-gaussian-kernel-using-fft#comment_896583

⋮## Direct link to this comment

https://ww2.mathworks.cn/matlabcentral/answers/547593-convolution-with-gaussian-kernel-using-fft#comment_896583

## Direct link to this comment

https://ww2.mathworks.cn/matlabcentral/answers/547593-convolution-with-gaussian-kernel-using-fft#comment_896613

⋮## Direct link to this comment

https://ww2.mathworks.cn/matlabcentral/answers/547593-convolution-with-gaussian-kernel-using-fft#comment_896613

Sign in to comment.