MATLAB: FFT nomalization problem

fftfft2fftshiftfourierMATLAB

Values of my intensity profile are much bigger after I used fftshift(fft2()) and values should be equal. My question is : is there a normalization problem when I do the fft() of my field U2 and what to do to have the right result ?
Here is my code so far :
 clear all;
 close all;
 clc;
 %%constants declaration and parameter calculation.
 E0 = 10;                                  %maximum field amplitude.
 z = 6.8e-3;                               %distance between fiber and lens in m.
 w0 = 15e-6;                               %waist of beamlets in m.
 lambda = 0.98e-6;                         %wavelength in m.
 Nx = 3;                                   %number of fiber over x axis.
 Ny = 3;                                   %number of fiber over y axis.
 dx = .4e-3;                               %pitche over x axis in m.
 dy = .4e-3;                               %pitche over y axis in m. 
 zr = pi*w0^2/lambda;                      %Rayleigh range in m.
 k = 2*pi/lambda;                          %wavenumber in m-1.
 wz = w0*sqrt(1+(z^2)/(zr^2));             %w(z) in m.
 Rz = z*(1+(zr^2)/(z^2));                  %R(z) in m.
 phiz = atan(z/zr);                        %Gouy phase.
 R = 1.4*wz;                               %lens radius in m. 
 %%[X,Y] domain creation.
 Npts = 2000;
 Nspots = 9;
 res1 = 1e-6;
 [X1,Y1] = meshgrid((-Npts/2:Npts/2-1).*res1,(-Npts/2:Npts/2-1).*res1);
 res2 = 50e-6;
 [X2,Y2] = meshgrid((-Npts/2:Npts/2-1).*res2,(-Npts/2:Npts/2-1).*res2);
 %%aperture center vectors.
 X0 = [-dx 0 dx -dx 0 dx -dx 0 dx];
 Y0 = [-dy -dy -dy 0 0 0 dy dy dy];
 %%field/intensity calculation and plot.
 U1 = zeros(Npts,Npts,Nspots);
 U2 = zeros(Npts,Npts,Nspots);
 tmp = E0*(w0/wz);
 for j = 1:length(X0)
     U1(:,:,j) = tmp*exp(-((X1-X0(j)).^2+(Y1-Y0(j)).^2)/(wz.^2));
     U2(:,:,j) = tmp*exp(-((X2-X0(j)).^2+(Y2-Y0(j)).^2)/(wz.^2));
 end
 U1 = sum(U1,3);
 U2 = sum(U2,3);
 UCC = conj(U1);
 I = U1.*UCC;
 figure(1);
 imagesc(X1(1,:),Y1(:,1),I);
 title('Real intensity profile');
 xlabel('x in m');
 ylabel('y in m');
 %%FFT calculation and plot.
 Uf = fftshift(fft2(U2));
 UfCC = conj(Uf);
 If = Uf.*UfCC;
 figure(2);
 imagesc(X2(1,:),Y2(:,1),If);
 title('FFT intensity profile');
 xlabel('x in m');
 ylabel('y in m');

Best Answer

I rarely use fft2, but I know that with fft it is necessary to normalise it by dividing it by the length of the time-domain signal.
See if:
Uf = fftshift(fft2(U2))/numel(U2);
does what you want.
Related Question