convolution with gaussian kernel using fft
74 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
David Goodmanson
2020 年 6 月 14 日
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 件のコメント
Paul Safier
2021 年 4 月 12 日
編集済み: Paul Safier
2021 年 4 月 12 日
@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 日
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.
その他の回答 (1 件)
Paul Safier
2021 年 4 月 14 日
編集済み: Paul Safier
2021 年 4 月 14 日
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 件のコメント
David Goodmanson
2021 年 4 月 22 日
編集済み: David Goodmanson
2021 年 4 月 22 日
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 日
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?
参考
カテゴリ
Help Center および File Exchange で Parametric Spectral Estimation についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!