current community

your communities

more stack exchange communities

company blog
Stack Exchange Inbox Reputation and Badges
tour help
Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote 12 down vote favorite
4

I have an optical pulse in time domain: 

Exp[-t^2] Cos[50 t - Exp[-2 t^2] 8 π].

The figure of this formula is

enter image description here

I hope to calculate the Fourier Transform of this formula, which gives the spectral distribution of this pulse. The spectral shape should be something looks like this:

enter image description here

But the Fourier Transform is difficult for the structure of Cos[t + Exp[t^2] ]. Anyone can help me by calculating the Fourier Transform of this formula? Or numerical Fourier Transform is also OK. Thank you very much!

share|improve this question
    
You might try it with the FourierSeries Package? – Phab 21 hours ago
    
Just out of curiosity, where do you find such an exotic pulse? – Jason B 20 hours ago
    
@Jason B, Thank you so much for your good answer. This pulse is from a self phase modulation effect in optical fiber, very common phenomenon in nonlinear optics. – user14634 20 hours ago
    
I'm curious to see what the spectrogram would look like for this pulse, I'm trying to work out a numeric Wigner function for it right now. – Jason B 20 hours ago
1  
To see the spectrum of this pulse, you can refer to the book: Nonlinear Fiber Optics (5th edition)-Govind P. Agrawal, 2013 – user14634 20 hours ago

3 Answers 3

up vote 14 down vote accepted

It always takes me a while to remember the best way to do a numerical Fourier transform in Mathematica (and I can't begin to figure out how to do that one analytically). So I like to first do a simple pulse so I can figure it out. I know the Fourier transform of a Gaussian pulse is a Gaussian, so 

pulse[t_] := Exp[-t^2] Cos[50 t]

Now I set the timestep and number of sample points, which in turn gives me the frequency range

dt = 0.05;
num = 2^12;
df = 2 π/(num dt);
Print["Frequency Range = +/-" <> ToString[num/2 df]];

Frequency Range = +/-62.8319

Next create a timeseries list, upon which we will perform the numerical transform

timevalues = 
  RotateLeft[Table[t, {t, -dt num/2 + dt, num/2 dt, dt}], num/2 - 1];
timelist = pulse /@ timevalues;

Notice that the timeseries starts at 0, goes up to t=num dt/2 and then goes to negative values. Try commenting out the RotateLeft portion to see the phase it introduces to the result. We will have to RotateRight the resulting transform, but it comes out correct in the end. I define a function that Matlab users might be familiar with,

fftshift[flist_] := RotateRight[flist, num/2 - 1];
Grid[{{Plot[pulse[t], {t, -5, 5}, PlotPoints -> 400, 
    PlotLabel -> "E(t)"],
   ListLinePlot[Re@fftshift[Fourier[timelist]], 
    DataRange -> df {-num/2, num/2}, PlotLabel -> "E(ω)"]}}]

enter image description here

which is what we were expecting. So now we try it on the more complicated pulse,

pulse[t_] := Exp[-t^2] Cos[50 t - Exp[-2 t^2] 8 π];
timelist = pulse /@ timevalues;
Grid[{{Plot[pulse[t], {t, -5, 5}, PlotPoints -> 400, 
    PlotLabel -> "E(t)"],
   ListLinePlot[Re@fftshift[Fourier[timelist]], 
    DataRange -> df {-num/2, num/2}, PlotLabel -> "Re E(ω)"]}}]

enter image description here

That doesn't look right ,the spectrum doesn't go to zero at the outer edges. We need more bandwidth on our transform, which we can get by decreasing the timestep

dt = 0.025;
df = 2 π/(num dt);
timevalues := 
  RotateLeft[Table[t, {t, -dt num/2 + dt, num/2 dt, dt}], num/2 - 1];
timelist = pulse /@ timevalues;
ListLinePlot[Re@fftshift[Fourier[timelist]], DataRange -> df {-num/2, num/2}, PlotLabel -> "Re E(ω)"]}}]

enter image description here

Or, if you want the power spectrum,

ListLinePlot[Abs@fftshift[Fourier[timelist]], DataRange -> df {-num/2, num/2}, PlotLabel -> "Abs E(ω)"]

enter image description here

share|improve this answer
    
Thanks a lot for this good answer! – user14634 20 hours ago
    
It would be the best if we can find a analytical answer. – user14634 20 hours ago
    
For that, maybe simplify the integral as much as you can and post it on the math board. – Jason B 20 hours ago
    
For reference: fftshift[dat_?ArrayQ] := Part[dat, Sequence @@ (Composition[RotateRight[#, Ceiling[Length[#]/2]] &, Range] /@ Dimensions[dat])] – J. M.♦ 8 hours ago
    
@Jason B I post the question here: math.stackexchange.com/questions/1533314/… – user14634 4 hours ago
up vote 10 down vote

Here's how to find the Fourier transform numerically (for a bandlimited signal).

Define the function (signal) of interest:

x[t_] := Exp[-t^2] Cos[50 t - Exp[-2 t^2] 8 \[Pi]];

Define the observation interval (this is necessarily finite, I used the values in your plot):

{ti, tf} = {-2, 2};

Now, we sample the signal with an appropriate sampling period ts. This must be small enough to avoid aliasing effects. Since I did not know the signal bandwidth beforehand, I had to find an appropriate sampling period by trial and error.

ts = 10.^-1.6;

The following list represents the sampled signal values:

xs = x /@ Range[ti, tf, ts];

The above list can now be transformed to the (discrete-time) frequency domain by means of Fourier function. Although you can directly apply Fourier to any list, I prefer to pad the latter with zeros so as to have an effective list length of $2^k,k\in \mathbb{N}$, as for that size the discrete fourier transform is more efficiently implemented (and you also gain in frequency resolution).

nfft = 1024;
xf = Fourier[PadRight[xs, nfft], FourierParameters -> {1, -1}];

I used the FourierParameters typically used in signal processing application, any other would do.

The required plot is the following (there are some details on how the xf points correspond to (continuous-time) frequencies, which I will not describe)

 ListPlot[Transpose[{Range[-nfft/2, nfft/2 - 1]/(nfft ts),
 RotateLeft[Abs[ts*xf], nfft/2]}], Joined -> True, PlotRange -> All, 
 Axes -> False, Frame -> True, FrameLabel -> {"f (Hz)", "magnitude"}]

enter image description here

share|improve this answer
    
Thanks a lot for this good answer! – user14634 20 hours ago
up vote 3 down vote

This is not the exhaustive answer, but the first step to it, after which you can proceed yourself, if you like this numerical approach. Try this:

f[t_] := Exp[-t^2] Cos[50 t - Exp[-2 t^2] 8 \[Pi]];        
tab = Table[
       NIntegrate[f[t]*Cos[k*t], {t, 0, \[Infinity]}, 
        Method -> {"LevinRule",  
          Method -> {"GaussKronrodRule", "Points" -> 11}}], {k, -300, 300,
         1}];

    ListLinePlot[tab, PlotRange -> All]

yielding a part of the Fourier-transform:

enter image description here

Have fun!

share|improve this answer
    
Thanks a lot for this good answer. But it takes me 30 seconds to run the program. – user14634 20 hours ago
    
This program is short and very useful. Can we make it faster? – user14634 20 hours ago
    
@user14634 I have no sure recipe. I would play with the parameters of the Method, and also try other methods suitable for the highly oscillating functions. – Alexei Boulbitch 14 hours ago
    
@ Alexei Boulbitch, thanks a lot! – user14634 4 hours ago

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.