convolution with gaussian kernel using fft

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]

Why aren't you using conv2()?

[LM]

It is an application to the theory of fast fourier transform, that's why I'm using fft.

Answer 1

[David Goodmanson]

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])

[LM]

Thanks a lot, now it's working! I just don't really understand your [3]: why should k needs a factor 2pi and what du you mean by w?

[David Goodmanson]

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]

That made everything clearer, thanks!

[Paul Safier]

This was a very helpful answer! Thank you @David Goodmanson!

[Paul Safier]

@David Goodmanson: Does n have to be odd, i.e. do you have to have a coordinate point identically at the origin? I have a kernel that goes as:

1/sqrt(x.^2 + y.^2)

and thus has a problem at the origin... I have added a small number (epsilon) to prevent the division by zero, but that adds error that might be avoided if I didn't have the exact origin in the domain. Any suggestions?

[David Goodmanson]

Hi Paul,

I'm not quite getting it. In the notation used to describe K at the top of the original question, does this mean 1/|x-y| ?

[Paul Safier]

@David Goodmanson Sorry, my question may have been confusing. I used your suggestions for the problem above on a different problem that I'm working with. I am solving the equation shown in the image attached below. I am able to successfully solve it with the following expressions (plus other supporting lines of code that I'm ommitting). The goal is to solve for p (press) in the equation since I know w and the kernel. The x/y domain can be the same as that of the original post where the origin is at the center.

% Set up the kernel
xfix = x; xfix(xfix == 0) = eps; % This removes zero from x so I never divide by zero.
Kern = 1./sqrt((xfix.^2+y.^2)); % Kernel in regular space. Pole at origin!
KernF = fftshift(fft2(ifftshift(Kern))); % fourier transform of kernel
wf = fftshift(fft2(ifftshift(w)));       % fourier transform of w
press = Const*fftshift(ifft2(ifftshift(wf./KernF)));

I have to add eps to prevent the discontinuity at the origin. Unfortunately my final result for press seems to depend on the magnitude of eps (eps ~ 1e-10 ish)! Do you have any advice for these situations?

[David Goodmanson]

Hi Paul,

I think there are some possibilities but I don't know of a good example with an analyical solution to compare to. Do you know of one?

[David Goodmanson]

Hi Paul,

here is one case with an analytical solution, which could be used for checks. In the original notation where x and y are 2d vectors

Int e^(-a*|x|^2) / |x-y| d2x 
    = pi^(3/2)/sqrt(a)  e^(-(a/2)*|y|^2) I0((a/2)*|y|^2)

where I0 is the modified bessel function of order 0, besseli.

Also if you go to fourier transforms,

    Int e^(-i*k.(x-y))/|x-y| d2x = 2*pi/|k|           (dot product)

where the integral can be d2x or d2y, doesn't matter.

[Paul Safier]

@David Goodmanson : Not sure I'm understanding. How come the left hand side (lhs) and right hand sides (rhs) aren't the same? Isn't the rhs a convolution?

x = linspace(-2,2,1000); y = x;
a = 2;
arg1 = a/2*abs(y).*abs(y);
lhs = pi^(3/2)/sqrt(a)*exp(-a/2*abs(y).*abs(y)).*besseli(0,arg1);
ff = exp(-a*abs(x).*abs(x));
gg = 1./abs(x);
rhs = conv(ff,gg,'same');
plot(y,lhs,'k',x,rhs,'g')

[David Goodmanson]

Hi Paul, the lhs is a convolution involving the (2d) integration variable x, and fiixed (2d) variable y. The rhs is not a convolution but rather a function of y, the result of the convolution.

[Paul Safier]

@David Goodmanson : but plots of the lhs and rhs does not lie on top of each other. They ought to be equal and when plotted, be on top of each other, no?

[David Goodmanson]

yes they should, but it is a 2d convolution.

[Paul Safier]

@David Goodmanson : Do you see what's wrong in this code to compare the left/right sides via a plot? I appreciate your help and this discussion...

x = linspace(-2,2,100); y = x; [xx,yy] = meshgrid(x,y);
a = 2;
arg1 = a/2*abs(yy).*abs(yy);
lhs = pi^(3/2)/sqrt(a)*exp(-a/2*abs(yy).*abs(yy)).*besseli(0,arg1);
ff = exp(-a*abs(xx).*abs(xx));
gg = 1./abs(xx);
rhs = conv2(gg,ff,'same');
plot(yy(:,5),lhs(:,5),'*k',xx(5,:),rhs(5,:),'g')

[David Goodmanson]

Hi Paul,

Your equations are not quite there yet for a 2d situation. Even though I used it a couple of times myself, the notation used for K in the original question [each of x and y is actually a 2d vector with x and y components] is not the best. The code below uses r with components rx and ry, and s with components sx and sy (although s is not really used).

The variable m picks a line of projection and compares the convolution (denoted c) and analytic (denoted cA) results as you were doing above. Agreement with theory is pretty good, and is the least good for a line that goes right across the singularity in the kernel, which is m = 201 here. If you zoom in for 201 you will see that the agreement in the peak value between c and cA is decent but not perfect.

The convolution has no singularity, but kernel g does at rx=ry=0. That can't be allowed to stand. For infinity I substituted the particular quantity 4/delx, where delx is the array spacing. I have a vague idea about how that is justifiable but have not proved it.

There is a factor of delx^2 in the convolution because conv is a 2d sum and the analytic calculation cA is a 2d integral.

The analytical solution contains the quantity

e^(-(a/2)*|r|^2) *I0((a/2)*|r|^2)

which is in the form

exp(-z)*besseli(0,z).

Not coincidentally that is exactly the normalized bessel function

besseli(0,z,1) 

so the code uses that.

N = 200;
delx = 1;
arr = (-N:N)*delx;
rx1 = arr;
ry1 = arr;
sx1 = arr;
sy1 = arr;
[rx ry] = meshgrid(rx1,ry1);
[sx sy] = meshgrid(sx1,sy1);
a = .2;
f = exp(-a*(rx.^2  +ry.^2));
%g = exp(-a*(sx.^2 + sy.^2));
g = 1./sqrt(sx.^2 +sy.^2);
g(isinf(g)) = 4/delx;
c = conv2(f,g,'same')*delx^2;
cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1);
% show some pictures
figure(1)
imagesc(arr,arr,log10(f))
figure(2)
imagesc(arr,arr,log10(g))
figure(3)
imagesc(arr,arr,log10(c))
figure(4)
imagesc(arr,arr,log10(cA))
m =  201;
xm = rx1(m);
figure(5)
plot(arr,c(m,:),arr,cA(m,:));
grid on
% get a bit away from the origin
m =  180;
xm = rx1(m);
figure(6)
plot(arr,c(m,:),arr,cA(m,:));
grid on

[Paul Safier]

I'll have a look. Thanks.

[Paul Safier]

@David Goodmanson : In toying with the code you give here, I see what was wrong earlier (thank you!), but I'm still left with my original confusion: the fidelity of the convolution to the analytical solution depends on the value that you replace infinity with in the kernel. In this example, I found the minimum error occurs when you use 3.8539 in the numerator instead of 4. Why this number? Intuitively it would seem that the larger number you use, short of infinity, would be most accurate since the true kernel is infinite at the origin. That is not the case though. If you use, for example, 100 in the numerator, your error is huge. Is this value a "known-good" for this specific kernel? If one has uneven grid spacing in x/y directions, which would you use: delx or dely, or an average of sorts? Do you know of any theory for this?

Testing of numerator:

num1 = linspace(.1,5,100); error = zeros(size(num1));
for ij = 1:length(num1)
g = 1./sqrt(sx.^2 +sy.^2);
g(isinf(g)) = num1(ij)/delx;
% g(isinf(g)) = 4/delx;
c = conv2(f,g,'same')*delx^2;
cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1);
error(ij) = cA(201,201) - c(201,201);
end
plot(num1,error)
c = polyfit(num1,error,1);
disp(['Equation is y = ' num2str(c(1)) '*x + ' num2str(c(2))])

[David Goodmanson]

@Paul Safier

Although the function g does have a singularity, it is not all that singular and goes away with the integration done in convolution. The 2d convolution is

c(r) = f(r-s)*g(s) d2s
where
g(s) = 1/|s|

with the singularity at sx = sx = 0. If you go to polar coordinates in a small area just around the origin then

sx = |s| cos(phi) 
sy = |s| sin(phi) 
dA = d2s = |s| ds dphi   % area element

plug that in get, for just the small-area-element part of the integral

c(r,dA) = f(r-s)*(1/|s|)*|s| ds dphi
     = f(r-s) ds dphi                       (1)

so the area element cancels the singularity, and using some very large number to approximate an infinity is not required.

You can use conv2 with its usual area element of dsx dsy for all but the one point, and make an approximation for that one point. Since (a) conv2 uses the multiplicative factor of delx^2, and (b) the area element in (1) has the single factor of ds which can be approximated by delx, you need to divide by delx for the approximation to work with conv. So the approx is const/delx.

Experimentation shows that const = 3.85 is good, and you even took it to a couple more decimal places, so with a bit more verification [ e.g. what happens if the constant 'a' in function f is no longer .2? ] that seems like a good value to use. Initially I had the idea that the constant might be 4, which is too large, and the latest value I found from theory is 3.53, which is too small but at least not far off.

Answer 2

[Paul Safier]

@David Goodmanson:

Can you point me to where you found the 3.53 value?

All this discussion has led me to conclude that since in my problem of interest, the unknown is within the convolution integral, when I manipulate the terms, I may be able to just multiply the kernel (in Fourier space) with the other unknown and remove the singularity. I seem to be off by a scale values (at least), do you know where?

As a test example, I would like to recover f in this expression (the equation we've been using), since I know g and c:

c(r) = f(r-s)*g(s) d2s

My strategy would be to form this expression (in essence):

f = ifft(fft(c)/fft(g))

But, since the analytical form of the fft of g is known, I believe:

g = 1/sqrt(x.^2 + y.^2);
fft(g) = 1/sqrt(kF.^2 + qF.^2); % According to Mathematica...

Then the quotient fft(c)/fft(g) is div1 below.

Here's the code:

N = 401; % Make odd
nov2 = (N-1)/2; 
delx = .25;
a = 0.2;
arr = (-nov2:nov2)*delx;
kk = (-nov2:nov2)*(1/delx)*2*pi/(N);
rx1 = arr; ry1 = arr; sx1 = arr; sy1 = arr;
[rx, ry] = meshgrid(rx1,ry1);
[sx, sy] = meshgrid(sx1,sy1);
[kF,qF] = meshgrid(kk,kk); % FFT space

f = exp(-a*(rx.^2  +ry.^2));
cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1);
cAF = fftshift(fft2(ifftshift(cA))); % FFT of analytical solution (RHS)
div1 = cAF.*sqrt(kF.^2 +qF.^2); % Here's the key part. This is fft(c)/fft(g). If I do this product, I can avoid the singularity. This needs a scalar multiplying it though??
ftest2 = 1/(delx*delx)*fftshift(ifft2(ifftshift(div1)));

% Compare
m =  201;
figure(5)
plot(arr,f(m,:),'k',arr,ftest2(m,:),'*');
grid on

Pretty sure somewhere I need some scale factors to toggle back and forth from Fourier space... Does this seem sensible??

[David Goodmanson]

Hi Paul,

< I think you have mistakenly posted the stuff above as an answer (easy to do, from personal experience) so if you wish to repost it as yet another comment to the comments I will save this and do the same. >

Dividing in the frequency domain is a bit trickier than multiplying, so I would recommmend as a first step, finding f(k) and g(k), then transforming their product to see if you can come up with cA. Making things harder for finding f by transforming cA/g is the fact that f(x) becomes impossibly small for larger values of x. It is going to become impossible to get an accurate f(x) for the full range of x.

For f(k)*g(k), again there will be a g singularity to deal with at kx = ky = 0.

The Mathematica result you quoted for the fft, let's say it is correct, but it leaves questions about intended normalization. For the fourier transform of g(x), a previous post of mine says, essentially,

Int e^(-i*k.x)/|x| d2x = 2*pi/|k|           (dot product)  

Since fft is being used to approximate an integral, I think it is safer to write down the integral and adjust the fft by multiplying by delx and so forth.

3.53: tbd

[Paul Safier]

@David Goodmanson : Oops, yes, meant to be a comment--must have hit the wrong button. Shall I delete and start this over?

Sorry, I'm having trouble understanding the notation:

Int e^(-i*k.x)/|x| d2x = 2*pi/|k|           (dot product) 

Is this saying the g(k) is 2*pi/|k|? Then is k = sqrt(kF.^2 + qF.^2)? Is the integral a convolution? Why are k (Fourier space) and x (normal space) both on the lhs?

Since g(x) is a fraction, the division cA/g is like a multiplicaton, no? Do you think this product is doomed?:

div1 = cAF.*sqrt(kF.^2 +qF.^2);

[Paul Safier]

Oh, I see, that's just the Fourier Transform integral... Sorry.

Another thing, in my problem, the scaling is such that x and y should never be greater than 10 and most likely less than 1...

[David Goodmanson]

since deleting an answer deletes all the comments to it (at least that is how it used to work, I don't know if that has changed), it may be too late to transfer over. But if you want to copy your answer over to the end of the original comment thread, then I copy over my reply, then you copy over your reply to that, etc then you delete the 'answer' and its comments, that could work.

[Paul Safier]

@David Goodmanson : Here's a test of a couple of possible ways of solving the problem I'm interested in. Can you let me know your thoughts? In my problem the unknown is inside the convolution integrand and this code uses your suggested test equation (with all three functions being known) and tries to recover f. Method1 evaluates the kernel with its singularity and then fixes it as we've been discussing. I look at a few values of the fix's numerator to gauge the sensitivity. Method2 multiplies the two functions (cA and g) to avoid dealing with the singularity. You mentioned issues with this approach, but I'm concerned with weighing those with the issue of using the "perfect" numerator for the fix. Which method is most susceptable to failing in a problem without a known analytical solution?

Running this will generate a plot comparing Methods 1 and 2 with the exact value of f.

% Test of methods of recovering f from within the convolution integral of the expression:
% g * f = cA
% f = exp(-a*(rx.^2  +ry.^2)) and g = 1./sqrt(sx.^2 +sy.^2) and cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1)

N = 401; % Make odd
nov2 = (N-1)/2; 
delx = .25;
arr = (-nov2:nov2)*delx;
xx = (-nov2:nov2)*delx;
kk = (-nov2:nov2)*(1/delx)*2*pi/(N);
rx1 = arr; ry1 = arr; sx1 = arr; sy1 = arr;
[rx, ry] = meshgrid(rx1,ry1);
[sx, sy] = meshgrid(sx1,sy1);
[kF,qF] = meshgrid(kk,kk); % FFT space
a = .5;
f = exp(-a*(rx.^2  +ry.^2)); % The function we are trying to recover

g = 1./sqrt(sx.^2 +sy.^2); % The kernel that has a singularity
g(isinf(g)) = 3.85/delx; % Replace the singularity
gF = fftshift(fft2(ifftshift(g))); % FFT of kernel
g1 = 1./sqrt(sx.^2 +sy.^2); % The kernel that has a singularity
g1(isinf(g1)) = 3.4/delx; % Replace the singularity
g1F = fftshift(fft2(ifftshift(g1))); % FFT of kernel
g2 = 1./sqrt(sx.^2 +sy.^2); % The kernel that has a singularity
g2(isinf(g2)) = 4.3/delx; % Replace the singularity
g2F = fftshift(fft2(ifftshift(g2))); % FFT of kernel

cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1);
cAF = fftshift(fft2(ifftshift(cA))); % FFT of analytical solution (RHS)
div1 = cAF.*sqrt(kF.^2 +qF.^2)*delx*delx/(2*pi);
ftestAlt = 1/(delx*delx)*fftshift(ifft2(ifftshift(div1)));
ftest = 1/(delx*delx)*fftshift(ifft2(ifftshift(cAF./gF)));
ftest1 = 1/(delx*delx)*fftshift(ifft2(ifftshift(cAF./g1F)));
ftest2 = 1/(delx*delx)*fftshift(ifft2(ifftshift(cAF./g2F)));

m =  nov2 + 1;
figure(1)
plot(arr,f(m,:),'k',arr,ftest1(m,:),'g',arr,ftest(m,:),'r',arr,ftest2(m,:),'b',arr,ftestAlt(m,:),'*');
xlim([-4 4]);
legend('Exact','Method 1 (num = 3.4)','Method 1 (num = 3.85)','Method 1 (num = 4.3)','Method 2')
grid on

[David Goodmanson]

Hi Paul,

I have not looked at your code yet, but when you say cA*g I assume what is meant is really cA/g.

When I said that division has some issues, consider convolving f and g to get h, with no assumptions yet about f,g,h. When you multiply f(k)*g(k) to get h(k) and transform to get h(t), that tends to work pretty well. But when you know g and h and want f, then of course f(k) = h(k)/g(k). Often g goes to 0 fairly quickly for large k. To get a good f, most of the time you need h to go to 0 about as quickly as g does. But if h is experimental data or something, it does not necessarily go to 0 as fast as it should. Then h/g gets large for large k, which messes up f.

[David Goodmanson]

Hi Paul,

I ran your code, and things look basically all right. I added the following

figure(6)
semilogy(arr,f(m,:),arr,ftest2(m,:));
grid on
figure(7)
semilogy(arr,f(m,:),arr,ftest2(m,:)/100);
xlim([-5.5 5.5])
grid on

Fig 6 shows that aside from being off by a normalization factor of approximately 100, the two agree pretty well for -5 <= k <= 5. Fig 7 is a blowup of fig. 6 wiith a factor of 100 thrown in. The agreement is good for about two decades of function value.

Although I am not anxious to do it, you can get the correct normalization by writing out the actual fourier integrals and then putting factors of delx^2 or delk^2 into the fft2 or ifft2. Be aware that for a 1d N-point ifft, the ifft does the sum and divides by N (that means divide by N^2 for a 2d ifft), so you have to correct for that.

As you can see the division f(k) = h(k)/g(k) has some problems when k gets too large and cA(k) gets too small.

[Paul Safier]

Hi @David Goodmanson :

I am confused as to the location of a factor of ~100 error you mention. I compared the fft/ifft Matlab summation formulae with the continuous integrals from the link below. It seems to suggest to convert the expressions from the commented to the uncommented shown below. The multiplying constant in front turns out to be numerically equivalent in both cases. However, I do not understand why it needed a second 1/(2 pi). I added that, willy-nilly, when I saw that the result was 2pi off of the exact solution! :(

cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1);
cAF = fftshift(fft2(ifftshift(cA))); % FFT of analytical solution (RHS)
div1 = cAF.*sqrt(kF.^2 +qF.^2)*delx*delx/(2*pi);

% ftestAlt = 1/(delx*delx)*fftshift(ifft2(ifftshift(div1)));
% ftest = 1/(delx*delx)*fftshift(ifft2(ifftshift(cAF./gF)));
% ftest1 = 1/(delx*delx)*fftshift(ifft2(ifftshift(cAF./g1F)));
% ftest2 = 1/(delx*delx)*fftshift(ifft2(ifftshift(cAF./g2F)));
dk = kk(2)-kk(1);
ftestAlt = N*N*dk*dk/(2*pi*2*pi)*fftshift(ifft2(ifftshift(div1))); % Method 2
ftest = N*N*dk*dk/(2*pi*2*pi)*fftshift(ifft2(ifftshift(cAF./gF))); % Method 1
ftest1 = N*N*dk*dk/(2*pi*2*pi)*fftshift(ifft2(ifftshift(cAF./g1F))); % Method 1
ftest2 = N*N*dk*dk/(2*pi*2*pi)*fftshift(ifft2(ifftshift(cAF./g2F))); % Method 1

https://reference.wolfram.com/language/ref/InverseFourierTransform.html

Bottom line is: I must solve for f somehow--my problem requires it. After I get to the bottom of the scaling issue, the question becomes, do I solve it via method 1 (replace the singularity) or method 2 (do the mult/div before evaluating the kernel (and its singularity)). I want to think method 2 is better since it doesn't rely on picking the appropriate constant to replace infinity with (since it seems arbitrary at least until I find theory discussing it). However, in messing with delx, N and such in the code above, I see that method 2 is less accurate in many instances. I don't see why that would be the case since both methods involve a division of two transformed terms. Why do div1 and cAF./gF lead to different final results??

Thanks for all your help and patience in working through this.

[David Goodmanson]

Hi Paul

kind of like a puzzle, not bad to work on. ok for the convolution we have

Int f(r')g(r-r') d2r' = h(r)             % convolution

f(k) = (1/(2*pi)) Int f(r)e^(-ik.r) d2r    % transform pair
f(r) = (1/(2*pi)) Int f(k)e^( ik.r) d2k

g(k) = (1/(2*pi)) Int g(s)e^(-ik.s) d2s    % transform pair
g(r) = (1/(2*pi)) Int g(k)e^( ik.s) d2k

h(r) =               Int h(k)*e^( ik.r) d2k  % transform pair
h(k) =  (1/(2*pi))^2 Int h(r)*e^(-ik.r) d2r
where h(k) = f(k)*g(k)

so f(r) = (1/(2*pi))*Int (h(k)/g(k))*e^( ik.r) d2k

For a 2d transform pair, you can pick any normalization constant c1 you want going into k space, but for the normalization constant c2 coming back the other direction, their product has to equal c1*c2 = (1/(2*pi))^2. Here I arbitrarily picked symmetric transform pairs for f and g so you don't have to remember anything extra. Then for the way the normalization works out, and defining h(k) = f(k)*g(k), the h transform pair is not symmetric.

If you define

delk = (1/delx)*2*pi/(N);
with
kk =  (-nov2:nov2)*delk;

and insert the appropriate factors (1/(2*pi)), delx^2 and delk^2 and N into the fft and ifft expressions, you should get agreement between the new f(r) and the original one. Let me know how that goes. You don't literally have to fourier transform g(s) since there is the analytic expression

g(k) = (1/(2*pi)) Int (1/|s|) e^(-i*k.s) d2s = 1/|k|  

and you don't have to do anything about g(k) at the origin since in the process |k| appears in the numerator so a zero is all right.

[Paul Safier]

Hi @David Goodmanson :

After running through that exercise, I keep getting an offset with the two methods. Method 1 replaces the singularity with 3.85/delx and method 2 avoids the singularity by multiplying cA(k)*|k|. Method 1 seems about dead-on with the known analytical solution for f(r). Method 2 comes up a bit short regardless of the mesh size. Do you know why?

This code compares the two methods against the analytical solution for f(r):

N = 801; % Make odd
xymax = 10.25;
nov2 = (N-1)/2; 
delx = xymax/nov2;
arr = (-nov2:nov2)*delx;
xx = (-nov2:nov2)*delx;
kk = (-nov2:nov2)*(1/delx)*2*pi/(N); dk = kk(2)-kk(1);
rx1 = arr; ry1 = arr; sx1 = arr; sy1 = arr;
[rx, ry] = meshgrid(rx1,ry1);
[sx, sy] = meshgrid(sx1,sy1);
[kF,qF] = meshgrid(kk,kk); % FFT space
a = .5;
f = exp(-a*(rx.^2  +ry.^2)); % The function we are trying to recover.

cA = (pi^(3/2)/sqrt(a))*besseli(0,(a/2)*(rx.^2+ry.^2),1);
cAF = fftshift(fft2(ifftshift(cA)));

% Method 1
g = 1./sqrt(sx.^2 +sy.^2); % The kernel that has a singularity
g(isinf(g)) = 3.85/delx; % Replace the singularity
gF = fftshift(fft2(ifftshift(g))); % FFT of kernel
meth1 = N^2*dk^2/((2*pi)^2)*fftshift(ifft2(ifftshift(cAF./gF)));

% Method 2 (uses analytical solution for fft of kernel)
meth2 = delx^2*N^2*dk^2/((2*pi)^3)*fftshift(ifft2(ifftshift(cAF.*sqrt(kF.^2 +qF.^2)))); % Avoid dealing with singularity

% Compare via a plot
m =  nov2 + 1;
figure(1)
plot(arr,f(m,:),'*k',arr,meth1(m,:),'g',arr,meth2(m,:),'r');
xlim([-0.4 0.4]);
legend('Exact','Method 1','Method 2')
grid on
% The ratio of the two methods is alway off!
meth1(m,m)/meth2(m,m)

[Paul Safier]

Here's a way I compared the arguments of both method's inversion:

% Test two arguments of methods 1 and 2
figure(3)
tstA = cAF./gF; % Method 1
tst1A = cAF.*sqrt(kF.^2 +qF.^2)*delx^2/(2*pi); % Method 2
plot(kk,abs(tstA(m,:)),'k',kk,abs(tst1A(m,:)),'g')
xlim([kk(m-50) kk(m+50)])
legend('Method 1','Method 2')

Looks like the origin (that was the location of the singularity) is the big difference. Does the origin need to be "fixed" in method 2??

[David Goodmanson]

Hi Paul,

[1] myself, I don't particularly like the line

cAF = fftshift(fft2(ifftshift(cA)))

because I start to lose track of normalization. The method I used is

cAF = (1/(2*pi))^2*delx^2*fftshift(fft2(ifftshift(cA))); % fourier trans of analytical solution (RHS)
% analytical result for fourier transform of g(r):
 gF = (1/2*pi) *Int e^(-ikr) / |r| d2r = 1/|k|
cAFdivgF = cAF.*sqrt(kF.^2 +qF.^2);
% fourier transform of the ratio
% meth3  = (1/(2*pi))*delk^2*N^2*fftshift(ifft2(ifftshift(cAFdivgF)));

(which gives the same comparison, and peaks out at .99).

But certainly any normalization you are good with is up to you. And it is true that if your immediate attention is not on frequency domain normalization, then a lot of the constants associated with mimicking integrals cancel out. If you define

cAF4 = fftshift(fft2(ifftshift(cA)));    (1)
cAFdivgF4 = cAF4.*sqrt(kF.^2 +qF.^2);    % cA(k)/g(k)
then
meth4 = [(1/(2*pi))*delk^2*N^2]*[(1/(2*pi))^2*delx^2]*fftshift(ifft2(ifftshift(cAFdivgF4)))
and since delx*delk = 2*pi/N,
meth4 = (1/(2*pi))*fftshift(ifft2(ifftshift(cAFdivgF)))  (2)

so (1) and (2) get the job done with a minimum of constants.

[2] as you recall, when I suggested replacing the 0 in g, it was for the convolution

Int f(r-s)*g(s) d2s.

Since there is a 0 of g(s) = 1/|s| in the denominator, replacing that 0 with something else is a necessity. For notation, suppose g0(s) stands for g without replacing the 0, and ge(s) is g, only having replaced the 0 with something else. So ge is well behaved, whereas g0 has a zero in the denominator and 1/g0 has a zero in the numerator. Since g0(s) = 1/|s| and g0(k) = 1\|k|, let the same notation apply in k space as well.

For the original convolution, you found ge replaces the 0 with is 3.85/delx, experimentally. So far so good,

For the deconvolution, with cA(k)/g(k), since there is a zero in the numerator the integral does not diverge and is well behaved, so there is no compelling reason to change the zero to anything else. I was quite surprised to find out that ge(k) works better than g0(k). Good catch by you. I have no idea at all why that might be the case, especially since you are replacing the zero of g0(k) -- in the k domain -- with 3.85/delx, which involves a spacing in the spacial domain.

(info on 3.53 to follow)

[David Goodmanson]

The 3.53 came about in the 2d convolution

h(r) = Int f(r-s)*g(s) d2s = Int f(r-s)*1/|s| d2s

Assume a 2d lattice of points for sx, sy, each with spacing delx. Draw a checkerboard of squares of side delx, with each point in the center of a square. Each point will represent the function within its square. For the square centered at the origin, assume that the square is small enough that f(r-s) is essentially constant. Going to polar coordinates u and phi, the integral over the square is

Int 1/|u| udu dphi = Int du dphi

and the singularity is gone. Assume phi = 0 is the x axis. The distance from the origin to the side is delx/2, and denote this by b. If the square were an inscribed circle of radius b, then the integral would just be 2*pi*b. But the square is a bit larger. Draw a 45 degree line from the origin to the upper right corner to make a triangle. There are eight such triangles. For this triangle, any ray from the origin to an intersection with the right hand side has length such that u*cos(phi) = b, so u = b/cos(phi). The u integration is

Int {0,b/cos(phi)}  du

which is just b/cos(phi) and the remaining integral is

Int {0,pi/4} b/cos(phi) dphi

Taking b outside temporarily,

syms phi
int(1/cos(phi))
ans = log(1/cos(phi)) + log(sin(phi) + 1)
 
int(1/cos(phi),0,pi/4)
ans = log(2^(1/2) + 1)
 
8*double(int(1/cos(phi),0,pi/4))  % 8 triangles
ans = 7.0510

The final result, 7.0510*b, is a bit larger than 2*pi*b.

Since b = delx/2, the final result is 3.53*delx which gets changed to 3.53/delx for the reason discussed before.

[Paul Safier]

Hi @David Goodmanson

Thanks for this explanation. I'm slightly confused still though since it seems it should be more straightforward as the area of a circle inscribed by that square would be pi/4*delx*delx and the area of the square would be delx*delx--both are different from the result obtained. It seems like the integration ought to provide an area, but its dimension is that of length and thus will shrink/grow as dx not dx^2. Curious. I'll think more about this. I imagine this issue must come up quite a bit in numerical analysis but I have not been able to find and literature regarding it. Do you know of a text, paper, web link or something else?

The difference between the methods is definitely perplexing. In my real problem that I'm dealing with, the results of the two methods are significantly different and it's definitely bothersome. However, on the other hand, by running that last code snippet I sent (that compares the arguments of both methods) it's obvious they are different and thus should NOT provide the exact same answer. But, why the "approximate" method (that replaces the singularity) is the one that's correct is the mystery...

Perhaps there's another analytical test case to look at. Did you find this example we've been looking at in a text/paper or did you derive it via analysis or symbolic convolution, etc?