% % Set up an initial Gaussian pulse and find its Fourier transform with FFT % % Input parameters tau0 = 10.e-12; % Assume Rayleigh length of 400m and require maximum propagation distance of 1km. % z0 will be replaced by the actual formula after we learn the Gaussian pulse z0= 400.; maxL=1.e3; % Estimate max. pulse width tau_max = tau0 * sqrt(1.+(maxL./z0).^2); trange = 5*tau_max sample= 50*trange./tau0; sample = 2^(fix(log(sample)./log(2.))+1) %convert it to 2^n for FFT deltaT=2*trange./sample; % Gaussian pulse with width tau0 and amplitude A index=0:1:(sample-1); t_shifted = index*deltaT; t = t_shifted - trange; frange = 0.5/deltaT; deltaf = 1 / (2*trange); f_shifted = index*deltaf; f = f_shifted - frange; A = exp(-(t./tau0).^2); t1 = t*1.e12; t1_shifted= t_shifted*1.e12; subplot(2,2,1); plot(t1,A); xlabel('time (ps)'); ylabel(['tau0 = ', num2str(tau0), 's']); title(['Gaussian Pulse in time domain']); axis([-trange*1.e12 trange*1.e12 0 1]); grid on; A = fftshift(A); %Prepare A in shifted format subplot(2,2,2); plot(t1_shifted,A); xlabel('time (ps)'); title(['Gaussian Pulse in time domain (shifted)']); axis([0 2*trange*1.e12 0 1]); grid on; Af = fft(A); % Fourier transform and discrete Fourier transform differ by deltaT % We also normalized it by tau0*sqrt(pi) Actual_Af = Af.*deltaT./(tau0.*sqrt(pi)); subplot(2,2,3); % Scale to THz f1_shifted =f_shifted*1.e-12 f1 = f*1.e-12 plot(f1_shifted,abs(Actual_Af)); xlabel('frequency (THz)'); ylabel('Normalized Amplitude (* tau0 * sqrt(pi))'); title(['Gaussian Pulse in frequency domain (shifted)']); axis([0 2.*frange*1.e-12 0 1]); grid on; subplot(2,2,4); plot(f1,abs(fftshift(Actual_Af))); xlabel('frequency (THz)'); title(['Gaussian Pulse in frequency domain']); axis([-1.e-12*frange 1.e-12*frange 0 1]); grid on;