MATLAB: Convolution avec noyau gaussien utilisant fft

Hey,
Je ne suis vraiment pas un pro en Matlab, j’ai donc quelques difficultés avec la tâche suivante.
J’essaie d’implémenter la diffusion d’un cercle par convolution avec le noyau gaussien 2D. La convolution est entre le noyau gaussien et la fonctionu, ce qui permet de décrire le cercle en étant +1 à l’intérieur du cercle et -1 à l’extérieur. Le noyau gaussien est
. J’ai essayé de ne pas utiliser fftshift mais de faire le shift à la main. Je sais aussi que la transformée de Fourier du gaussien estavec des coefficients dépendant de la longueur de l’intervalle.
Mais avec mon code, il n’y a aucune diffusion et je ne comprends pas ce que j’ai fait de mal. J’apprécierais vraiment de l’aide.
J’ai le code suivant avec un maillage 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;
 

Meilleure réponse

  • Salut LM
    Le code ci-dessous reprend votre approche mais modifie certains détails. [1] Les tailles de tableau sont impaires x impaires puisque vous obtenez ainsi la plus grande symétrie. [2] Comme vous pouvez le voir, le code utilise beaucoup de fftshift et ifftshift. Pour les tableaux impairs, fftshift n’est pas son propre inverse et ifftshift n’est pas non plus son propre inverse. Mais ce sont des inverses les uns des autres. Tout ce dont vous devez vous souvenir, c’est que fftshift prend k = 0 ou x = 0 au centre du tableau, et ifftshift les ramène tous les deux à l’origine afin que fft et ifft fonctionnent correctement. [3] le tableau k a besoin d’un facteur de 2pi par rapport à ce que vous avez, puisque k = 2pi (1 / lambda), analogue à w = 2pi f.
    [4] J’ai utilisé une fonction qui est 0 en dehors du cercle au lieu de -1. L’utilisation de -1 donne un énorme pic à k = 0 dans le domaine k, mais l’utilisation de 0 donne un bien meilleur tracé de ce qui se passe dans le domaine k. Vous pouvez passer à -1 au point indiqué dans le code et le résultat est bon pour ce cas. [5] Vous n’avez pas à vous soucier de la normalisation du noyau dans l’espace k. C’est parce que la fonction est exp (-k. ^ 2 * t). Vous savez que la fonction devrait = 1 lorsque t = 0, car il ne devrait y avoir aucun effet. La constante devant doit donc être 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])