convolution with gaussian kernel using fft

Question

LM 2020-6-13

1
链接

此问题的直接链接

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

评论： Paul Safier 2021-4-22

采纳的回答： David Goodmanson

在 MATLAB Online 中打开

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;

2 个评论
显示无隐藏无

Image Analyst 2020-6-13

Why aren't you using conv2()?

LM 2020-6-13

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2020-6-14

3
链接

此回答的直接链接

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

在 MATLAB Online 中打开

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

18 个评论
显示 16更早的评论隐藏 16更早的评论

Paul Safier 2021-4-1

编辑：Paul Safier 2021-4-1

在 MATLAB Online 中打开

@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 2021-4-9

编辑：David Goodmanson 2021-4-9

在 MATLAB Online 中打开

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 2021-4-12

编辑：Paul Safier 2021-4-12

在 MATLAB Online 中打开

@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 2021-4-13

编辑：David Goodmanson 2021-4-13

在 MATLAB Online 中打开

@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 2021-4-14

0
链接

此回答的直接链接

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

编辑：Paul Safier 2021-4-14

在 MATLAB Online 中打开

@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??

14 个评论
显示 12更早的评论隐藏 12更早的评论

David Goodmanson 2021-4-14

编辑：David Goodmanson 2021-4-14

在 MATLAB Online 中打开

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 2021-4-15

编辑：Paul Safier 2021-4-15

在 MATLAB Online 中打开

@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

Paul Safier 2021-4-19

在 MATLAB Online 中打开

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 2021-4-19

编辑：David Goodmanson 2021-4-20

在 MATLAB Online 中打开

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 2021-4-20

在 MATLAB Online 中打开

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)

David Goodmanson 2021-4-21

编辑：David Goodmanson 2021-4-21

在 MATLAB Online 中打开

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 2021-4-22

编辑：David Goodmanson 2021-4-22

在 MATLAB Online 中打开

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 2021-4-22

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?

请先登录，再进行评论。

convolution with gaussian kernel using fft

2 个评论
显示无隐藏无

采纳的回答

18 个评论
显示 16更早的评论隐藏 16更早的评论

更多回答（1 个）

14 个评论
显示 12更早的评论隐藏 12更早的评论

另请参阅

类别

标签

Community Treasure Hunt

convolution with gaussian kernel using fft

2 个评论 显示 无隐藏 无

采纳的回答

18 个评论 显示 16更早的评论隐藏 16更早的评论

更多回答（1 个）

14 个评论 显示 12更早的评论隐藏 12更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

18 个评论
显示 16更早的评论隐藏 16更早的评论

14 个评论
显示 12更早的评论隐藏 12更早的评论