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

Gives me

enter image description here

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...

$\endgroup$
2

2 Answers 2

16
$\begingroup$

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]
$\endgroup$
5
  • $\begingroup$ Great answer. How do I get the right values on the x-axis? I understand that that the DataRange option would be appropriate. And that first I have to select only the first half of the graph. But then I also need to know what the maximum frequency is... which may be simple in this case but what if I didn't know what function I did the Fourier transform on? $\endgroup$
    – C. E.
    Jan 20, 2013 at 8:27
  • 2
    $\begingroup$ @Pickett The maximum frequency will be SamplingFrequency/2 and if you have an arbitrary function where you don't know if the end points line up, you typically use a window function to bring the end-points to zero. It sounds like you're having trouble in the signal processing part than Mathematica. If so, I would suggest Signal Processing for some of the conceptual questions on DSP :) $\endgroup$
    – rm -rf
    Jan 20, 2013 at 16:05
  • $\begingroup$ Perhaps it's a bit of both, yes. Anyway a lot of great help here. Thank you a lot. $\endgroup$
    – C. E.
    Jan 20, 2013 at 16:20
  • 1
    $\begingroup$ @rm-rf Indeed, but you have to take care that in order to get the maximum frequency at the end of the scale you have to drop the second half of the list. This also prevents the second peak from showing up. DataRange can be used to get the axis labeling right. Set it to {0, SamplingFrequency/2}. $\endgroup$ Jan 20, 2013 at 18:48
  • $\begingroup$ @Sjoerd C. de Vries: Thanks for the $DataRange trick. Now how to adjust the amplitude? amplitude comes with a multiplier of 20. why it is so? $\endgroup$ Sep 19, 2013 at 18:17
10
$\begingroup$

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.

$\endgroup$
1
  • 5
    $\begingroup$ Could you add corresponding images? $\endgroup$
    – ybeltukov
    Sep 20, 2013 at 13:20

Your Answer

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

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