0
$\begingroup$

I am trying to compute the Fourier transform of a list, say an array of the form {{t1, y[t1]},.....{tn, y[tn]}}; apply some filters in the spectral components, and then back transform in time domain.

I have written my own routine to do Discete Fourier Transform of lists, as follow

FourierAmpl[data_] := 
 Module[{ydata, xdata, \[CapitalDelta]freqdata, xfreqdata, datatime, 
   datafourier},
  ydata = 
   Join[Transpose[data][[2]], 
    Table[Transpose[data][[2]][[-1]], {ii, 1, zerotabbing}]]; 
  xdata = Join[Transpose[data][[1]], 
    Table[Transpose[data][[1]][[-1]] + \[CapitalDelta]x ii, {ii, 1, 
      zerotabbing}]]; \[CapitalDelta]freqdata = 1/(
   xdata[[Length[xdata]]] - xdata[[1]]); 
  xfreqdata = 
   Table[\[CapitalDelta]freqdata ii, {ii, 0, Length[xdata] - 1}]; 
  datafourier = 
   Drop[Transpose[{xfreqdata, Abs[Fourier[ydata - Mean[ydata]]]/
      1}], -Round[0.5 Length[xfreqdata]]];
  datafourier
  ]

which directly computes the Fourier Transform of a 2D list, provided that these quantities are defined zerotabbing and \[CapitalDelta]x, which respectively indicate a number of extra zeros that can be added after the signal, to improve the spectral resolution; while \[CapitalDelta]x the spacing in the time signal (dt).

I have also written a function to compute also the argument of the FFT:

renormalize[args_] := 
 Module[{pairs, diffs, j, len = Length[args], corr = 0}, 
  pairs = Partition[args, 2, 1];
  diffs = Map[#[[1]] - #[[2]] &, pairs];
  PrependTo[diffs, 0];
  diffs = 2*Pi*Sign[Chop[diffs, Pi]];
  Table[corr += diffs[[j]];
   corr + args[[j]], {j, 1, len}]]

FourierPhase[data_] := 
 Module[{ydata, xdata, \[CapitalDelta]freqdata, xfreqdata, argdata, 
   phasedata},
  ydata = 
   Join[Transpose[data][[2]], 
    Table[Transpose[data][[2]][[-1]], {ii, 1, zerotabbing}]]; 
  xdata = Join[Transpose[data][[1]], 
    Table[Transpose[data][[1]][[-1]] + \[CapitalDelta]x ii, {ii, 1, 
      zerotabbing}]]; \[CapitalDelta]freqdata = 1/(
   xdata[[Length[xdata]]] - xdata[[1]]); 
  xfreqdata = 
   Table[\[CapitalDelta]freqdata ii, {ii, 0, Length[xdata] - 1}];
  argdata = renormalize[Arg[Fourier[ydata - Mean[ydata]]]];
  phasedata = 
   Drop[Transpose[{xfreqdata, argdata}], -Round[
      0.5 Length[xfreqdata]]];
  phasedata
  ]

The code seems to work in properly computing the FFT. For example, for a very simple function such as

Test = Table [{tt, Sin[2 \[Pi] tt]}, {tt, 0, 6, 0.05}]; \[CapitalDelta]t = Test\[Transpose][[1]][[2]] - Test\[Transpose][[1]][[1]]; zerotabbing = 0;

Test function

The function

ListLinePlot[FourierAmpl[Test], PlotRange -> All]

returns

FFT of Test

But now if I would like to recompute back the original signal applying twice FourierAmpl the resulting function is far from the original function:

ListLinePlot[FourierAmpl[FourierAmpl[Test]], PlotRange -> All]

enter image description here

Where is my mistake? I fail in finding it. Thanks for your help!

Improve this question
$\endgroup$
3
  • $\begingroup$ Did you consider the built-in function? If you do this: test = Table[{tt, Sin[2 [Pi] tt]}, {tt, 0, 6, 0.05}]; testDFT = Fourier[test[[;; , -1]], FourierParameters -> {1, -1}]; Then this: ListLinePlot[Thread[{test[[;; , 1]], Abs[testDFT]^2}], PlotRange -> All] and ListLinePlot[ Thread[{test[[;; , 1]], Re[InverseFourier[testDFT, FourierParameters -> {1, -1}]]}]] you get back where you started. I realize this isn't pointing out the mistake in your code. I do see that the value returned from the built-in is different than yours. $\endgroup$
    – Mark R
    Commented Jun 21, 2019 at 21:28
  • $\begingroup$ If the forward FFT has a plus sign in the exponent Exp[I w t] then the inverse needs a minus sign Exp[-I w t]. There may also be a scaling issue, but this may be the source of your problem. $\endgroup$
    – bill s
    Commented Jun 22, 2019 at 2:03
  • $\begingroup$ Thanks a lot @MarkR. I also saw that the built in function work. My problem with that approach, is that I have lots of dataset with different spacing in dt and different length. Also I will need to apply filters in the frequency domain. This will make it very difficult to compare different signals if I am using the built in ones $\endgroup$
    – Matteo S
    Commented Jun 22, 2019 at 8:02

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.