- 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.
- 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.

# The difference between performing the fft2 and the convolution

4 views (last 30 days)

Show older comments

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

##### 0 Comments

### Answers (2)

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.

Hope this helps!

##### 0 Comments

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

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!