How can I get a Fourier-limited Gaussian pulse

27 views (last 30 days)
Hi all, I have generated a femtosecond Gauss pulse and fft it to get the frequency spectrum. But when I calculate the time bandwidth product, it is never a Fourier limited pulse, e.g., the time duration (FWHM) is 20fs, and the bandwidth is around 0.2eV(read from the graph), and 20fs*0.2eV is never close to the theoretical limit 1.825 eV*fs. In principle, this Gaussian pulse be Fourier limited. I don't know if I make some mistakes here. I hope that someone can help me with this. Thanks! Below is my code:
t=[-100:0.01:100];
A=0.1519;
T0=20; %FWHM of the pulse
w=2*pi*300/800; %frequency centered at 800nm
y=A.*exp(-4*log(2)*(t/T0).^2).*cos(w.*t); %Gauss pulse in terms of FWHM expression
figure;plot(t,y);
Y=fftshift(fft(y));
Y_mag=abs(Y);
X=(-(length(Y_mag)-1)/2:(length(Y_mag)-1)/2)*1239.84/300/length(Y_mag)/(t(2)-t(1)); %frequency axis
figure;plot(X,Y_mag);

Accepted Answer

David Goodmanson
David Goodmanson on 7 Sep 2018
Hi Andi,
The issue here is that for a voltage pulse the fourier transform is done on the linear voltage function y(t) to produce Y(f), but the FWHM condition is on the square of y(t) and Y(f). I modified your code to produce a gaussian pulse y(t) as the square root of your y(t) [which in this context is actually y(t)^2]. Then keep the cos(wt) carrier, do the fft to get Y(E), and plot y^2 and Y^2 for FWHM. The height of the peaks don't matter for FWHM, so I made no attempt to figure out the constants in front.
For the Y(E)^2 plot, the FWHM of each peak is about .09 eV (it's a bit undersampled), and 20fs times that is 1.8, so it works.
I added figure numbers since I'm not a fan of new figures popping up every time you rerun the code.
t=-100:0.01:100;
A=0.1519;
T0=20; %FWHM of the pulse
w=2*pi*300/800; %frequency centered at 800nm
y=sqrt(exp(-4*log(2)*(t/T0).^2)).*cos(w.*t); %Gauss pulse in terms of FWHM expression
y2 = y.^2;
figure(1)
plot(t,y2);
Y=fftshift(fft(y));
Y2_mag=abs(Y.^2);
X=(-(length(Y_mag)-1)/2:(length(Y_mag)-1)/2)*1239.84/300/length(Y_mag)/(t(2)-t(1)); %frequency axis
figure(2)
plot(X,Y2_mag,'-o');
xlim([1.4 1.7])
ylim([0 3e6])
  3 Comments
David Goodmanson
David Goodmanson on 9 Sep 2018
It's partly due to historical reasons. Suppose you have a decaying exponential signal
y(t) = e^(iw0t)e^-(g/2)t t >= 0 g stands for gamma
(a complex signal is easier than using, for example) cos(w0t) ).
Without worrying about constants in front, its fourier transform is
Y(w) = 1/(w-w0-ig/2)
Its absolute value is
1/sqrt((w-w0)^2 +g^2/4)
Mathematically it is usually easier to deal with abs^2 of something rather than abs of something. Abs^2 is
1/((w-w0)^2 +g^2/4)
This peaks out at w = w0, and has half its value at w = w0 +- (g/2). The FWHM = g, so you can theoretically read g off of squared frequency peaks. You could get the same result by dropping down by 1/sqrt(2) on abs(Y(w)), but this is usually not the way it is done.
Same argument in time domain, although FWHM in that case does not usually have as direct an interpretation.
LL
LL on 10 Sep 2018
Got it, thank you very much, David!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!