Phase in discrete Fourier transformation

Question

I asked a question some days ago and it was very well answer. However, I have a new one. I don't understand why the phase of the discrete fourier transformation is not correct. It has the correct answer, but it also has other points that are not correct.

Here, the correct analytical answer:

FourierTransform[Exp[-t]*UnitStep[t], t, w, 
 FourierParameters -> {1, -1}]
Plot[Abs[-(I/(-I + 2*Pi*f))], {f, -5, 5}, PlotRange -> All]
Plot[Arg[-(I/(-I + 2*Pi*f))], {f, -5, 5}, PlotRange -> All]

Here, the numerical results I obtained:

signal[i_] := Exp[-(i*samplingPeriode)]*UnitStep[(i*samplingPeriode)]

initialTime = -10;
finalTime = 10;
signalDuration = finalTime - initialTime;

bandwidth = 50;
nyquistRate = 2*bandwidth;
nyquistSamplingPeriode = 1/nyquistRate;

samplingFrequency = 2*nyquistRate;
samplingPeriode = 1/samplingFrequency;
numberOfSamples = signalDuration/samplingPeriode;

yValuesSampledSignal = 
  Table[signal[i], {i, initialTime/samplingPeriode, 
    finalTime/samplingPeriode,samplingPeriode}];
xValuesSampledSignal = 
  Table[i*samplingPeriode, {i, initialTime/samplingPeriode, 
    finalTime/samplingPeriode,samplingPeriode}];
sampledSignal = 
  Partition[Riffle[xValuesSampledSignal, yValuesSampledSignal], 2];

yValuesDiscreteFourierTransform = 
  samplingPeriode*
   Fourier[yValuesSampledSignal, FourierParameters -> {1, -1}];
rotatedYValuesDiscreteFourierTransform = 
  RotateRight[yValuesDiscreteFourierTransform, numberOfSamples/2];
xValuesDiscreteFourierTransform = 
  Table[i*(1/signalDuration), {i, -(numberOfSamples/2), 
    numberOfSamples/2}];
discreteFourierTransform = 
  Partition[
   Riffle[xValuesDiscreteFourierTransform, 
    yValuesDiscreteFourierTransform], 2];

absYValuesDiscreteFourierTransform = 
  Abs[rotatedYValuesDiscreteFourierTransform];
absDiscreteFourierTransform = 
  Partition[
   Riffle[xValuesDiscreteFourierTransform, 
    absYValuesDiscreteFourierTransform], 2];

argYValuesDiscreteFourierTransform = 
  Arg[rotatedYValuesDiscreteFourierTransform];
argDiscreteFourierTransform = 
  Partition[
   Riffle[xValuesDiscreteFourierTransform, 
    argYValuesDiscreteFourierTransform], 2];

yValuesRecontructedSignal = (1/samplingPeriode)*
   InverseFourier[yValuesDiscreteFourierTransform, 
    FourierParameters -> {1, -1}];
xValuesRecontructedSignal = xValuesSampledSignal;
reconstructedSignal = 
  Partition[
   Riffle[xValuesRecontructedSignal, yValuesRecontructedSignal], 2];

ListLinePlot[sampledSignal, PlotRange -> All]
ListPlot[Take[
  absDiscreteFourierTransform, {numberOfSamples/2 - 100, 
   numberOfSamples/2 + 100}], Filling -> Axis, PlotRange -> All]
ListPlot[Take[
  argDiscreteFourierTransform, {numberOfSamples/2 - 100, 
   numberOfSamples/2 + 100}], Filling -> Axis, PlotRange -> All]
ListLinePlot[reconstructedSignal, PlotRange -> All] 

Answer 1

A few points:

• You need to work on giving the minimal working example that reproduces your problem if you can.
• To that end, your example here uses a timestep of 0.000025, leading to a sampled signal with eight hundred thousand elements. This makes it take an unreasonably long amount of time to work on and plot. You could make an example with only four thousand points that would give exactly the same plots.
• Whenever possible, use numeric quantities instead of integers. In your code, you can do this by changing initialTime, finalTime, and bandwidth to be -10.0, 10.0, and 50.0, respectively. This cuts down the execution time by half on my machine.
• It is confusing that you always use Partition[ Riffle[ a, b], 2] when Transpose[ {a,b}] is much shorter - but this is just a matter of style.
• When you define discreteFourierTransform, you use yValuesDiscreteFourierTransform, but when you define absYValuesDiscreteFourierTransform you use rotatedYValuesDiscreteFourierTransform

The actual issue at hand is that you need to rotate your time series before applying the DFT, so that the zero time is first. The xValuesSampledSignal should be define as (letting timerange/2 and -timerange/2 be the max and min time values, and dt be the timestep)

{0, dt, 2 dt, .....timerange/2,-timerange/2 + dt,-timerange/2 + 2 dt,...-2 dt, -dt}

To show this, here is a reworking of your code (to shorten it), first using a linear time list:

(* Define the time range and the number of sample points, this 
determines the frequency range and resolution. In order for the 
fftshift function to work properly, the number of sample points 
should be even.  *)
trange = 20.0;
num = 2^12;
dt = trange/num;
df = 1/(num dt);


fftshift[flist_] := RotateRight[flist, (Length@flist)/2 - 1];

signal[t_] := Exp[-t] UnitStep[t]

(* Now define the sampled signal in the time domain *)

xValuesSampledSignal := Table[t, {t, -dt num/2 + dt, num/2 dt, dt}];
yValuesSampledSignal = signal /@ xValuesSampledSignal;

(* Then define the spectrum, and plot the signal along with the amplitude and phase of the spectrum *)
yValuesDiscreteFourierTransform = 
  fftshift[Fourier[yValuesSampledSignal, 
    FourierParameters -> {1, -1}]];
Row[{ListLinePlot[
   Transpose[{xValuesSampledSignal, yValuesSampledSignal}], 
   PlotLabel -> "Signal"],
  ListLinePlot[Abs[dt yValuesDiscreteFourierTransform], 
   PlotLabel -> "Abs spectrum", DataRange -> df {-num/2 + 1, num/2}, 
   PlotRange -> {{-5, 5}, All}], 
  ListLinePlot[Arg[dt yValuesDiscreteFourierTransform], 
   PlotLabel -> "Arg spectrum", PlotRange -> {{-5, 5}, All}, 
   DataRange -> df {-num/2 + 1, num/2}]}]

This is all wrong, so now we try it with the shifted time values.

xValuesSampledSignal := 
  RotateLeft[Table[t, {t, -dt num/2 + dt, num/2 dt, dt}], num/2 - 1];
yValuesSampledSignal = signal /@ xValuesSampledSignal;
yValuesDiscreteFourierTransform = 
  fftshift[Fourier[yValuesSampledSignal, 
    FourierParameters -> {1, -1}]];

Row[{ListLinePlot[Abs[dt yValuesDiscreteFourierTransform], 
   PlotLabel -> "Abs spectrum", DataRange -> df {-num/2 + 1, num/2}, 
   PlotRange -> {{-5, 5}, All}], 
  ListLinePlot[Arg[dt yValuesDiscreteFourierTransform], 
   PlotLabel -> "Arg spectrum", PlotRange -> {{-5, 5}, All}, 
   DataRange -> df {-num/2 + 1, num/2}]}]

Notice that there is a slight slope in the Arg plot, due to the fact that the values at the endpoints need to match (that is, the phase goes to the same value at $\pm \omega_{\mathrm{max}})$. I used 4096 data points in these plots, if you change this to 800000 as in your example, this effect is not visible.