The difference between performing the fft2 and the convolution

4 views (last 30 days)
Hello,
I am trying to get the wavevector distribution of the field EP (ES*ES) via the fft2 simply or the convolution method. The distinct different results puzzled me. The ES represents a Gaussian beam with normalized power, and EP is the electrical field distribution of another beam (ES*ES). Naturally, we can get the wavevector distribution by performing the fft2, and the obtained wavevector distribution is consistent with that of the therotical case. Accroding to the convolution theorem, the wavevector amplitude of the beam EP can be described by the convolution of the fft(ES). In the following code, a loop is applied to implement the convolution. Each component of the EP's wavevectors (EP_K(fy,fx)) is a sum of the many terms of the product with the fft(ES) (ES_K(fy0,fx0)*ES_K(fy-fy0,fx-fx0)). In fact, there is a term, formed as ' sqrt(a-fx0^2-fy0^2)-sqrt(b-(fx-fx0)^2-(fy-fy0)^2) ', which is ignored in the attached code. This term changes for the different product of the fft(ES) although the (EP_K(fy,fx)) keep unchanged, which forced me to take a circular computation. However, the results from the convolution calculation is inconsistent with the fft2 computation.
And any suggestions to reduce the time spent on calculating this?
I appreciate your help very much. Thanks in advance!
Below is the code.
clear;
alpha=5;
Nx=128;
Ny=128;
D0=1e-3; %%%% the beam diameter
omega_0=D0/2; %%%% the beam radius
delta_x=2*alpha*omega_0/Nx;
delta_y=2*alpha*omega_0/Ny;
delta_xy=delta_x*delta_y;
x=(-Nx/2:1:Nx/2-1)*delta_x;
y=(-Ny/2:1:Ny/2-1)*delta_y;
[X,Y]=meshgrid(x,y); %%%% real space
kxmax=1/delta_x;
kymax=1/delta_y;
kx=(-(Nx)/2:1:(Nx)/2-1)*kxmax/(Nx);
ky=(-(Ny)/2:1:(Ny)/2-1)*kymax/(Ny);
delta_kx=kxmax/(Nx);
delta_ky=kymax/(Ny);
[KX,KY]=meshgrid(kx,ky); %%%% space frequency domain
ES=sqrt(2/pi)/omega_0*exp(-(X.^2+Y.^2)/omega_0^2); %%%% the electrical field of a gaussian beam with normalized power
ES_K=fftshift(fft2(fftshift(ES)));
ES_K=ES_K*delta_xy;
E1T=ES.*conj(ES);
E1T=sum(sum(E1T))*delta_xy; %%%%%%% the power in real space
E1t=ES_K.*conj(ES_K);
E1t=sum(sum(E1t))*delta_kx*delta_ky; %%%%%%% the power in wavevector space
ES_Kd=omega_0*sqrt(2*pi)*exp(-omega_0^2*4*pi^2*(KX.*KX+KY.*KY)/4); %%%% the deduced wavevector distribution
E1td=ES_Kd.*conj(ES_Kd);
E1td=sum(sum(E1td))*delta_kx*delta_ky;
%%%%%%%%%The calculation of the wavevector distribution via the fft method or the convolution method
EP=2/pi/omega_0^2*exp(-2*(X.*X+Y.*Y)./omega_0^2); %%%%% EP=ES*ES
EP_Kfft=fftshift(fft2(fftshift(EP)))*delta_xy;
EP_Kd=exp(-omega_0^2*4*pi^2*(KX.*KX+KY.*KY)/8); %%%% the deduced wavevector distribution
kx0=delta_kx;
ky0=delta_ky;
EP_Kconv=zeros(Ny,Nx);
for count_lefty=1:Ny
for count_leftx=1:Nx
P_plinshi=0;
for count_midy=1:Ny
for count_midx=1:Nx
count_rigy=count_lefty-count_midy+Ny/2+1;
count_rigx=count_leftx-count_midx+Nx/2+1;
if count_rigy<=Ny && count_rigy>=1 && count_rigx<=Nx && count_rigx>=1
P_plsls=ES_K(count_midy,count_midx)*ES_K(count_rigy,count_rigx); %%%%accroding to the convolution theorem
P_plinshi=P_plinshi+P_plsls;
end
end
end
EP_Kconv(count_lefty,count_leftx)=P_plinshi;
end
end

Answers (2)

UDAYA PEDDIRAJU
UDAYA PEDDIRAJU on 24 Nov 2023
Hi Chen,
I understand the discrepancies you faced between the results obtained from “fft2” and convolution methods to acquire the wavevector distribution of the field “EP (ES*ES)”. Also, I understand that you’re trying to optimize the computational time. I suggest the following approaches to resolve the issues.
  1. Vectorization: Utilize MATLAB's vectorization capabilities to optimize the loop-based convolution calculation for improved computational efficiency. you can refer to https://in.mathworks.com/matlabcentral/answers/2032009-what-is-vectorization-why-use-vectorization-instead-of-loops.
  2. Utilize MATLAB's Built-in Functions: Leverage MATLAB's built-in functions for convolution and Fourier transforms, such as “conv2” and “fft2”, which are optimized for efficiency. You can refer to https://www.mathworks.com/help/matlab/ref/fft2.html and https://www.mathworks.com/help/matlab/ref/conv2.html for more details.
Hope this helps!

Bruno Luong
Bruno Luong on 24 Nov 2023
Use this function, it will perform the convolution using fft method rightly
The syntax is more or less compatible with MATLAB conv2

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!