4
$\begingroup$

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]

enter image description here

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] 

enter image description here

$\endgroup$
2
  • $\begingroup$ Try yValuesSampledSignal = Table[signal[i], {i, 1, signalDuration/samplingPeriode}]; xValuesSampledSignal = Table[i*samplingPeriode, {i, 1, signalDuration/samplingPeriode}]; $\endgroup$
    – andre314
    Dec 26, 2015 at 23:29
  • $\begingroup$ This seems to solve it, but it eliminates the part before time=0 of the function, and I want to preserve this part of the function. $\endgroup$
    – Gabriel
    Dec 27, 2015 at 0:17

1 Answer 1

4
$\begingroup$

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}]}]

enter image description here

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}]}]

enter image description here

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.

$\endgroup$
3
  • $\begingroup$ Great answer. I think it will be nice if you provide some explanation on why you need to rotate the time series before DFT. $\endgroup$
    – BlacKow
    Dec 29, 2015 at 0:51
  • $\begingroup$ Thank you Jason! As BlacKow, why so complicated? Why so many rotations? And why the phase goes to the same value at ±ωmax? It makes it weird. and with more complicated signals it is very confusing. $\endgroup$
    – Gabriel
    Dec 30, 2015 at 1:40
  • $\begingroup$ @BlacKow and Juan, I'm not really an expert on DFTs, I just use them and know that I have to do it like this or it won't work. In fact, I often have to remind myself exactly the right way to do it by taking the DFT of a simple Gaussian whose transform I know, in order to get the parameters all correct. This is a good place to start for understanding why you need to set up the input data in this palindrome-like structure. $\endgroup$
    – Jason B.
    Dec 30, 2015 at 9:27

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