Fourier optics
HW 4

define functions for normalized FFTs in the t and ω domains

FFT[fLst_,npts_,dt_]:=(Sqrt[npts (dt^2)])RotateLeft[Fourier[RotateRight[fLst,npts/2]],npts/2];
InverseFFT[FLst_,npts_,dt_]:=(Sqrt[(1/(npts (dt^2)))])RotateLeft[InverseFourier[RotateRight[FLst,npts/2]],npts/2];

Problem 1

Using the FFT and convolve.nb program as a guide, create a pulse function, put it on a grid, and calculate and plot the final pulse shape for pure 2nd, 3rd and 4th order spectral phase. Do this for a Gaussian and a rectangular pulses of FWHM of 10fs and 100fs and a center wavelength of 800nm. 
Symbolically: start with (Subscript[E, in])(t), transform to (Subscript[E, in])(ω), multiply by Exp[ i ϕ(ω) ], then inverse FT to (Subscript[E, out])(t), then plot ((Abs[ (Subscript[E, out])(t) ])^2)

The Taylor expansion of spectral phase is defined as: 
ϕ(ω)=(Subscript[ϕ, 0])+(Subscript[ϕ, 1])(ω-(Subscript[ω, 0]))+(1/2)(((Subscript[ϕ, 2])(ω-(Subscript[ω, 0])))^2)+(1/3)(((Subscript[ϕ, 3])(ω-(Subscript[ω, 0])))^3)+(1/4)(((Subscript[ϕ, 4])(ω-(Subscript[ω, 0])))^4)+...,
so specifically, you’ll be plotting the output pulse shape for different values of (Subscript[ϕ, 2]), (Subscript[ϕ, 3]), (Subscript[ϕ, 4]), with (Subscript[ϕ, 0])=(Subscript[ϕ, 1])=0.
It is conventional to represent the (Subscript[ϕ, n])'s in units of (fs^n). 

So you should make 6 plots for each of the two pulse shapes: 
τ = 10fs: (Subscript[ϕ, 2]), (Subscript[ϕ, 3]), (Subscript[ϕ, 4])
τ = 100fs: (Subscript[ϕ, 2]), (Subscript[ϕ, 3]), (Subscript[ϕ, 4])
When plotting the temporal pulse shapes, determine the amount of each type of spectral phase that approximately triples the duration of the pulse. For each plot, make it clear what value of phase you are using, and plot the x-axis in units of fs. 
To do this, set up a grid for the time points of your array, and do a ListPlot[ Transpose[ {timeLst / fs, data} }, Joined  -> True], potentially using PlotRange to zoom in to get a better display.  
Keep an eye on whether the pulse broadens so much that the edges get close to the boundaries of the grid. If it does, you can increase the number of points or increase the dt spacing.

problem 2:

A signal f(t)=(sinc^2)((Subscript[ω, 1]) t), where (Subscript[ω, 1])=100, is sampled by a comb function, to produce the list (Subscript[f, s])(t)=f(t)comb((Subscript[ω, s])t), where (Subscript[ω, s]) is the sampling frequency (in rad/s).

part a:

Find the analytic transform, F(ω)=FFT{f(t)} and make plots of f(t) and F(ω). Use the internally defined function Sinc[ ].

part b:

Find the minimum sampling rate (Subscript[ω, s])=(Subscript[ω, Ny]) (Nyquist rate) and the corresponding sampling interval (Subscript[t, Ny])=2π/(Subscript[ω, Ny]) that will permit an exact recovery of f(t) from (Subscript[f, s])(t); i.e. such that g(t)=f(t) within a multiplicative constant.

part c:

Generate a list (Subscript[f, s])(t) (using the Table[ ] function in Mathematica) with a sampling period (Subscript[t, s])=2π/(Subscript[ω, Ny]) for 32 points.  

Use the format of ListPlot[ ] shown below to plot (Subscript[f, s])(t). These options for ListPlot[ ] allow you to see dots for the points as well as connecting the lines.

ListPlot[fInLst,Mesh->All,Joined->True,MeshStyle->PointSize[Medium],PlotRange->All]

Take the FFT (using the FFT[ ] function defined in fft and convolve.nb), and plot the spectrum. 

Repeat the above for 64 points with the same sampling interval (Subscript[t, s])=2π/(Subscript[ω, Ny]). Comment on what the effect is in both domains when you increase npts while maintaining dt.

part d:

As we showed in the class notes, the analytic calculation of the FT of the sampled function results in the convolution of the spectrum with another comb function. In terms of equations: FT{ (Subscript[f, s])(t) } = FT{ f(t) comb[t/(Subscript[t, s])] } = (Subscript[F, s])(ω)=(1/(2π))F(ω)⊗[(Subscript[ω, s])comb((ω/(Subscript[ω, s])))]. 
Aside: to obtain the same amplitude scaling for this result, so that the peak value of (Subscript[F, s])(ω) equals the peak value of F(ω), I find that we need to divide by dω=((2π)/(N dt))=((Subscript[ω, s])/N), so that
(Subscript[F, s])(ω)=(1/(2π dω))F(ω)⊗[(Subscript[ω, s])comb((ω/(Subscript[ω, s])))]=(N/(2π))F(ω)⊗[comb((ω/(Subscript[ω, s])))]

Make a plot of this analytic (Subscript[F, s])[ω] over a range from -2(Subscript[ω, s]) to 2(Subscript[ω, s]). To represent the comb[ω/(Subscript[ω, s])] function, you can use the following definition,

comb[ω_,ωs_]:=Sum[DiracDelta[ω-n ωs],{n,-3,3}];

then you can use the Mathematica function Convolve[ ] to do the analytic convolution. 

How does the result of the FFT numerical transform relate to the analytic result?

part e:

Repeat the plots made in parts c and d above, but for under-sampling (Subscript[t, s])=(3/2)((2π)/(Subscript[ω, Ny])),  and over-sampling, (Subscript[t, s])=(1/2)((2π)/(Subscript[ω, Ny])).

part f:

We now want to see what the recovered (continuous) function g(t) is like. In the frequency domain, applying a filter function H(ω)=rect[ω/(Subscript[ω, s])] selects the central copy of the spectrum. In the time domain, we can use the invFT of H(ω) as the interpolating function to reconstruct the signal. 
Find the filter impulse response h(t), and make a plot of it. 

Next, find the reconstructed function g(t) by convolving h(t) with (Subscript[f, s])(t). Since h(t) is a continuous function, you can’t use ListConvolve[ ]. Instead, create an interpolated function by summing up the products: g(t)=(∑_(-N/2)^(N/2))(Subscript[f, s])(n dt) h(t-n dt). 
Plot the reconstructed function g(t) for the case where sampling is at the Nyquist limit.

part g:

If you return to the case where (Subscript[ω, s])=(Subscript[ω, Ny]), note that the actual values of the sampled points depends on the phasing of the sampling. To illustrate this, shift f(t) by (Subscript[t, s])/2, then sample it to get (Subscript[f, s])(t). In the end, the recovered function g(t) should be the same. Show this is true by plotting the new g(t).