Plotting the frequency spectrum of a data series using Fourier

Question

testData = Table[N@Sin[500 x], {x, 0, 100}];
ListLinePlot[Abs[Fourier[testData]], PlotRange -> Full]

Gives me

Which I do not expect because the Fourier Transform is FourierTransform[Sin[500 x], x, f],

I Sqrt[Pi/2] DiracDelta[-500+f]-I Sqrt[Pi/2] DiracDelta[500+f]

I'm not saying that Mathematica's Fourier function is somehow faulty. But could someone please explain why there isn't a peak around 500? This is the frequency of the signal after all...

Answer 1

As mentioned by @Rahul, you have not sampled your sine wave often enough and have introduced artifacts due to aliasing. The frequency of Sin[500 x]=Sin[2 Pi f x] is $f=500/(2\pi)$, which is about 80 Hz. At least two samples per cycle are required to avoid aliasing, hence the default $x$ interval of 1 in {x,0,100} must be reduced to less than about $1/(2*80)=0.006$.

In addition, the discrete fast Fourier transform assumes periodicity. Hence, care must be taken to match endpoints precisely. An interval without an exact integral multiple of the sine wavelengths will return blurred Dirac delta functions.

I set the sampling interval to $(1/f)/4$, which is small enough to avoid aliasing. There are an integral number of wavelengths, 100*(2 Pi/500)/4, and the endpoints are matched.

testData = Table[N@Sin[500 x], 
              {x, 0, 100*(2 Pi/500)/4 - (2 Pi/500)/4, (2 Pi/500)/4}];
ListLinePlot[Abs[Fourier[testData]], PlotRange->Full]

Answer 2

This is some code I wrote long before features like DataRange were added to MMA. It is part of a larger package I wrote a long time ago and it is far (really far) from elegant. Perhaps the OP can find them useful and with some rewriting could even make them look better on newer version of MMA.

Purpose of these snippets is to automate the scaling of the time and frequency axes based on one of the parameters supplied as an option to the AddXXXRange procedures. Starting time and starting frequency are understood to be both zero.

Options[Domains] ={
          DeltaT -> None, TotalT -> None, DeltaF -> None, TotalF -> None,
      Toffset -> 0, Centered -> False,OnlyParameters -> False};

Domains[n_, opts___] :=
  Module[
    {tos, Ts, Ttot, Fs, Ftot,
      centered, param, timedomain, frequencydomain},

  Ts = DeltaT /. {opts} /. Options[Domains];
  Ttot = TotalT /. {opts} /. Options[Domains];
  Fs = DeltaF /. {opts} /. Options[Domains];
  Ftot = TotalF /. {opts} /. Options[Domains];
  tos = Toffset /. {opts} /. Options[Domains];
  centered = Centered /. {opts} /. Options[Domains];
  param = OnlyParameters /. {opts} /. Options[Domains];

  Switch[
      First[Flatten[Position[N[{Ts, Ttot, Fs, Ftot, tos}], _?NumberQ]]],
      1, Ttot = n Ts; Fs = 1/Ttot; Ftot = 1/Ts,
      2, Ts = Ttot/n; Fs = 1/Ttot; Ftot = 1/Ts,
      3, Ts = 1/(n Fs); Ttot = 1/Fs; Ftot = n Fs,
      4, Ts = 1/Ftot; Ttot = n Ts; Fs = Ftot/n,
      _, Ts = 1; Ttot = n; Fs = 1/n; Ftot = 1
  ];

  If[param,
      Return[{Ts, Ttot, Fs, Ftot}],
      timedomain = tos + (Range[n] - 1)Ts;
      frequencydomain =
          If[centered,
              (Range[n] - Ceiling[n/2])Fs,
              (Range[n] - 1)Fs
          ];
      Return[{timedomain, frequencydomain}]
  ]
]

Options[AddSignalRange] = {Toffset -> 0};

AddSignalRange[data_List, opts___] :=
  Transpose[{Domains[Length[data], opts][[1]], data}]

Options[AddSpectrumRange] = {
      Centered -> False,
      HalfSpectrum -> False
      };

AddSpectrumRange[data_List, opts___] :=
  Module[
      {n, lista, half, centered},

      half = HalfSpectrum /. {opts} /. Options[AddSpectrumRange];
      If[half,
          centered = False,
          centered = Centered /. {opts} /. Options[AddSpectrumRange];
      ];

      n = Length[data];
      lista = If[centered,

              Transpose[
                          {Domains[n, Centered -> True, opts][[2]],
                              RotateRight[data, Ceiling[n/2] - 1]}
              ],

              Transpose[
                          {Domains[n, Centered -> False, opts][[2]],
                              data}
              ]   
      ];
      If[half, (*then centered is forced to False*)
                      Take[lista, {1, Ceiling[(n + 1)/2]}],
                      lista (*can be centered or not *)
  ]
]

Now, let's take an example signal with components at 4 and 12 Hz, and let's sample it with a sampling frequency of 30 Hz. If we want 150 sample points, we will span a time interval going from 0 to 149 1/Fs. Fs being the sampling frequency.

n = 150; Fs = 30;
data = Table[.7Sin[2Pi 4 t] - .4Cos[2Pi 12 t] + .2(Random[] - .5),
        {t, 0, (n - 1)/Fs, 1/Fs}];

When we plot the naked data, we have n on the x-axes

ListLinePlot[data, PlotRange -> Full]

We can add the correct time scale with AddSignalRange, by specifying the sampling frequency used to gather that data

signal = AddSignalRange[data, TotalF -> Fs];
ListLinePlot[signal, PlotRange -> Full]

Same goes with the FFT of the sample data. We can perform our operation on the naked data

mag = Abs[Fourier[data]];

and then we apply the correct frequency range using AddSpectumRange. We can also specify we want an unfolded spectrum with negative and positive frequencies

spectrum = AddSpectrumRange[mag, TotalF -> Fs, Centered -> True];
ListLinePlot[spectrum, PlotRange -> Full]

or, just get the positive half

spectrum = AddSpectrumRange[mag, TotalF -> Fs, HalfSpectrum -> True];
ListLinePlot[spectrum, PlotRange -> Full]

The peaks are at 4 and 12 hertz.