18
$\begingroup$

I looked at all the other questions related to mine before posting this, and they didn't solve my problem. I have a large data set which can be downloaded from here. I'm using the following code to plot the DFT:

data = ToExpression @ Import["http://ge.tt/api/1/files/3dvZMs61/0/blob?download"];
ListPlot[Abs[Fourier[data]], Joined -> True, AxesOrigin -> {0, 0}]

enter image description here

As far as I can tell, there is no pronounced single or countable number of frequencies in this data set. However, if I just plot the data, I get:

ListPlot[data, AxesOrigin -> {0, 0}]

enter image description here

Here it is evident that there is a subtle periodic character, whose frequency is what I want. Why am I not getting that frequency as a peak in the first plot?

Improve this question
$\endgroup$
2
  • $\begingroup$ You can try removing the biggest outlier. $\endgroup$
    – b.gates.you.know.what
    Dec 3, 2013 at 19:17
  • 2
    $\begingroup$ @aman_abhishk: When I first started programming in Mathematica, I ran into this same "problem" and was very confused for a while! The problem is the plot range. Unless told otherwise, Mathematica automatically determines a vertical plot range which shows a lot of squiggles on screen, to put it nonrigorously. As a result, your y-axis is way too zoomed in, so things look weird. Try running your code again, but use ListPlot[Abs[Fourier[data]], Joined -> True, PlotRange -> All], and you'll see what's actually going on. $\endgroup$
    – DumpsterDoofus
    Dec 26, 2013 at 15:32

2 Answers 2

Reset to default
19
$\begingroup$

update

Just to clean things up a bit, we can use the discussion here to make a couple functions that help extract the frequency data from this dataset. I define two functions findPeriod and reconstruct:

Clear[findPeriod];
findPeriod[data_, threshold_] := 
 Module[{fs, s1, s = {}, i, a0f, af, pf, pos, fr, frpos, fdata, 
   fdatac, n, per},
  n = Length[data];
  fs = Fourier[data];
  s1  = Drop[fs, -Floor[Length[fs]/2]];
  For[i = 1, i < Length[s1], i++, 
   If[Abs[fs][[i + 1]] > threshold, AppendTo[s, i + 1]]];
  a0f = Abs[fs[[1]]]/Sqrt[n];
  af = 2/Sqrt[n] Abs[fs][[s]];
  pf = Arg[fs][[s]];
  {a0f, Transpose[{s, af, pf}]}]

Clear[reconstruct];
reconstruct[data_, fp_] := Module[{n},
  n = Length[data];
  Show[
   ListLinePlot[data, PlotStyle -> Black],
   Plot[fp[[1]] + 
     Sum[fp[[2, j, 2]] Cos[
        2 Pi (fp[[2, j, 1]] - 1)/n t - fp[[2, j, 3]]], {j, 1, 
       Length[fp[[2]]]}], {t, 0, n}, PlotStyle -> Red]
   ]]

The first can be used to extract all frequencies above a certain threshold:

TableForm[Sort[findPeriod[data, 0.4][[2]], #1[[2]] > #2[[2]] &], 
 TableHeadings -> {None, {"Freq", "Amp", "Phase"}}]

Mathematica graphics

and the second takes this information to reconstruct the data:

enter image description here

old version

I had to drop the first and last points of your TXT file, so there may be some slight differences. Process should be the same, however.

f = Abs[Fourier[data]];
n = Length[f];
(* n = 3737 *)
(* Drop the first point, which may be an outlier, find the max *)
pos = Position[f, Max[f[[2 ;;]]]][[1, 1]];
(* pos = 40 *)
(* Follow the help file for Fourier *)
fr = Abs[Fourier[data Exp[2 Pi I (pos - 2) N[Range[0, n - 1]]/n], 
    FourierParameters -> {0, 2/n}]];
frpos = Position[fr, Max[fr]][[1, 1]];
(* frpos = 2790 *)
(* Get the period *)
per = N[n/(pos - 2 + 2 (frpos - 1)/n)];
(* per = 94.6252 *)
(* Plot the extracted frequency with the data *)
pdata = Table[Mean[data] + 0.02 Sin[2 \[Pi] x/per], {x, n}] ;
ListLinePlot[{pdata, data}]

Mathematica graphics

Improve this answer
$\endgroup$
3
  • $\begingroup$ I don't know how or if the offset (which I used as Mean[data] and the amplitude (which I fudged) can be obtained. $\endgroup$
    – bobthechemist
    Dec 3, 2013 at 20:31
  • $\begingroup$ +1 for the update. I was doing just that (albeit along a different route). This is the second time in this post that I arrive late :-), the first being a comment (later deleted) with the same content as Simon Wood's post. How I hate being old and slow! $\endgroup$
    – Peltio
    Dec 4, 2013 at 22:07
  • $\begingroup$ @Peltio I'm usually on that end, so thanks for letting me beat someone for once. $\endgroup$
    – bobthechemist
    Dec 4, 2013 at 23:27
17
$\begingroup$

Looking at your plotted data you can see about 40 cycles of the dominant frequency, this tells you that the peak will appear somewhere around the 40th element of the DFT. That's in the region where your plot of the DFT is clipped, so it's no wonder you can't see the peak.

Looking at the relevant part of the DFT you can see the peak quite clearly:

ListLinePlot[Abs[Fourier[data][[2 ;; 100]]], PlotRange -> All, DataRange -> {2, 100}]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ You pointed out a very interesting point that if there is N cycles of the dominant frequency, this tells you that the peak will appear somewhere around the Nth element of the DFT. I would like to ask a question that "if we consider a wider initial range, i.e., we will have more cycles, the peak will have a larger value, right?". So when we compute the Fourier spectrum of a time-varying signal, the frequency of the peak we find will depend on the duration we consider but not on the frequency of the signal? So how can we find the frequency of the signal? Thank you very much! Tiao $\endgroup$
    – user42410
    Aug 17, 2016 at 16:49
  • $\begingroup$ @Tiao, with a signal of longer duration the DFT has greater frequency resolution. The peak appears at the same frequency, but that frequency appears at a different position in the output list. You can probably find a detailed explanation on the mathematics or signal processing sites, it's not a Mathematica specific thing. $\endgroup$
    – Simon Woods
    Aug 17, 2016 at 20:54

Your Answer

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

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