# How can I use Mathematica to numerically compute a Wigner spectrogram of an optical pulse?

This question was inspired by this question where it is necessary to numerically compute the Fourier transform of a Gaussian optical pulse with a Gaussian chirp function.

$$E(t)=e^{-t^2} \cos(50 t - e^{-2 t^2} 8 π)$$.

I was thinking that the Fourier transform, while a complete picture of the pulse, wasn't as useful as I would want it to be. I think a spectrogram, showing the time and frequency information of the pulse, would be helpful as well. The only spectrogram I am familiar with is the Wigner function (useful reading here as well)

$$W_x(t,f) = \int_{-\infty}^\infty E\left(t+\frac{\tau}{2}\right) E^*\left(t-\frac{\tau}{2}\right) e^{-i 2\pi\tau\,f} \, d\tau$$

What is the best way to implement this in Mathematica?

My first attempt was to do it directly,

ω0 = 50;
pulse[t_] := Exp[-t^2] Exp[-I (ω0 t - 8 π Exp[-2 t^2] )];
pulsecc[t_] := Exp[-t^2] Exp[I (ω0 t - 8 π Exp[-2 t^2] )];
wignerfunc[t_, w_] :=
NIntegrate[(pulse[t + τ/2] pulsecc[
t - τ/
2]) Exp[-I w τ], {τ, -∞, ∞}];
Plot[wignerfunc[0, w], {w, 30, 80}]

But I get convergence errors,

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in τ near {τ} = {-0.925724}. NIntegrate obtained 2.42861*10^-17+8.97719*10^-17 I and 3.6945011190775015`*^-15 for the integral and error estimates. >>

Is there an option to give NIntegrate that makes this converge?

This paper describes a discrete Wigner transform, that I imagine would take a 2D array (t along one dimension, τ along the other) and spit out the spectrogram, but signal processing always seems opaque to me.

• If you want a spectrogram, you could also consider ContinuousWaveletTransform. – Rahul Nov 18 '15 at 11:47
• @Rahul, I have no experience with that function, I'm interested to look at it. Also, I wonder how Spectrogram would compare to the result below. – Jason B. Feb 4 '16 at 19:28

This is the best I can come up with, I'm very interested to see if anyone else has a better solution. The idea here is to just run through values of $t$, and do a DFT on

$$E(t+\frac{\tau}{2}) E ^*(t-\frac{\tau}{2})$$

So I set up the time/frequency resolution for my DFT, using a dt value I know gives a broad enough spectrum,

dt = 0.025;
num = 2^14;
df = 2 π/(num dt);
timevalues =
RotateLeft[Table[t, {t, -dt num/2 + dt, num/2 dt, dt}], num/2 - 1];
fftshift[flist_] := RotateRight[flist, num/2 - 1];
Print["Frequency Range = +/-" <> ToString[num/2 df]];

Then I define the kernel of the transformation function

ω0 = 50;
pulse[t_] := Exp[-t^2] Exp[-I (ω0 t - 8 π Exp[-2 t^2] )];
pulsecc[t_] := Exp[-t^2] Exp[I (ω0 t - 8 π Exp[-2 t^2] )];
wignerkernelfunc[t_, τ_] := (pulse[t + τ/2] pulsecc[t - τ/2])

Then I build up the table , which is , doing a DFT at every time step,

Monitor[spectrogram =
Chop @ Table[timelist = wignerkernelfunc[t, #] & /@ timevalues;
fftshift[Fourier[timelist]], {t, -2, 2, .01}];, t]

And then I plot the distribution function, zooming in on the relevant region of the spectrum,

parulacolorlist = Get["http://pastebin.com/raw.php?i=HjYmaRRq"];
ParulaCM = Blend[Apply[RGBColor, parulacolorlist, {1}], #1] &;
ListDensityPlot[
spectrogram[[All, Round[15/df + num/2] ;; Round[85/df + num/2]]],
DataRange -> {{15, 85}, {-2, 2}}, PlotRange -> All,
ColorFunction -> ParulaCM, PlotLegends -> Automatic,
FrameLabel -> {Style["\[Omega]", 20], Style["t", 20]}]

You can plot 1D slices of this for single $t$ values,

ListLinePlot[
spectrogram[[{201 - 20, 201, 201 + 20},
Round[15/df + num/2] ;; Round[85/df + num/2]]],
DataRange -> {15, 85},
PlotLegends -> {"t = -0.2", "t = 0", "t = 0.2"},
FrameLabel -> {"ω", "amplitude"}]

or for constant ω, and here I vertically offset the plots for clarity,

ω1 = 40;
ω2 = 50;
ω3 = 60;
ListLinePlot[
Evaluate[Transpose[
spectrogram[[;; ,
Round[#/df +
num/2] & /@ {ω1, ω2, ω3}]]] + {0, .7,
1.4}], AspectRatio -> 1,
PlotLegends ->
Evaluate["ω = " <>
IntegerString[#] & /@ {ω1, ω2, ω3}],
FrameLabel -> {"t", "amplitude"}, DataRange -> {-2, 2}]

Edit: In the absence of an actual discrete Wigner transform, as described in the signal processing paper above, I think this is a good solution. I was able to produce basically the same plot using NIntegrate, via

wignerfunc[t_, w_] :=
NIntegrate[(pulse[t + τ/2] pulsecc[t - τ/2]) Exp[
I w τ], {τ, -∞, ∞}];
spectrogram2 = Table[Quiet@wignerfunc[t, w], {w, 15, 85, .25}, {t, -2, 2, .01}]

But it took over 10 times as long, and I had to use Quiet to stop getting convergence errors. I would be happy to hear what Method to use for NIntegrate to handle this better though.

I would also be interested to see what other forms of spectrogram would look like for a pulse like this. I don't quite understand what ContinuousWaveletTransform does. Another route would be to use a sliding window Fourier transform.

• Your approach appears fine to me. I am not an expert, but I believe the "checkerboard" features are typical in Fourier-based time-frequency analysis due to the inherent limit in time-frequency resolution. – Stelios Nov 17 '15 at 16:02
• @Jason B, This Wigner function is very helpful to understand the full picture in both time and frequency domains. I tried to run your program, but after the counts are finished, there is no figure output. At the same time, the computer is almost down (system halted). What is the reason? – user14634 Nov 20 '15 at 9:50
• @user14634, do you get a numerical output from spectrogram[[1]]? Also, did you do the top part of the code or the bottom part, which uses NIntegrate? – Jason B. Nov 20 '15 at 10:04
• @Jason B, I used your program from dt = 0.025; num = 2^14; df = 2 π/(num dt); ...... to ListDensityPlot[spectrogram......]. The output is Frequency Range = +/-125.664 – user14634 Nov 25 '15 at 7:42
• @Jason B In[1]: Get["pastebin.com/raw.php?i=HjYmaRRq"] Out[1]=\$Failed warnings: Get::noopen: "Cannot open \!(\"pastebin.com/raw.php?i=HjYmaRRq\"). " – user14634 Nov 25 '15 at 7:49