The difference between performing the fft2 and the convolution

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

Answer 1

[UDAYA PEDDIRAJU]

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!

Answer 2

[Bruno Luong]

Use this function, it will perform the convolution using fft method rightly

https://www.mathworks.com/matlabcentral/fileexchange/24504-fft-based-convolution

The syntax is more or less compatible with MATLAB conv2