MathGroup Archive 1999

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: searchforperiod

  • To: mathgroup at smc.vnet.net
  • Subject: [mg19771] Re: searchforperiod
  • From: Paul Abbott <paul at physics.uwa.edu.au>
  • Date: Fri, 17 Sep 1999 01:36:34 -0400
  • Organization: University of Western Australia
  • References: <7rnkao$elp@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

Andre Hautot wrote:

> I have computed the time evolution of a certain quantity, say rmod, (the
> details of the physical problem which leads to them are unimportant).
> The results are contained in a list like this :
> Table[{t[i],rmod[i]},{i,0,5000}]
> To fix the ideas here is the beginning of a typical list :
> {{0,1.},{2.02484567313,1.30384048104},{5.44775639416,2.34594834824},
> {11.921842635,4.27850743029},{23.2431295354,7.13760085727},
> {41.2272127052,10.9087937072},{67.6554674978,15.562780576},
> {104.248787003,21.0622014538},{152.649762838,27.3628881652},
> {214.407985645,34.4143233727},{290.966976416,42.1600354445},
> {383.652320768,50.53803291},{493.660902853,59.4812881541},
> {622.051243359,68.918265984},{769.734982324,78.7734911707},...}
> The rmod-values increase during 35 steps and they decrease during the 35
> next steps, returning near the initial value of one after 70 steps, and
> so on. Similarly, the time spacing increases during 35 steps, decreases
> during the 35 next steps and so on.
> The coordinates are known with arbitrary high precision (50 figures for
> example or more if you need)
> The graph of rmod versus t is obtained, as usual, by
> ListPlot[Table[{t[i],rmod[i]},{i,0,5000}]]
> It seems to be periodic. How can verify such a conjecture and obtain a
> high precision value for the period?

One approach is to use the Lomb normalized periodogram
(see Numerical Recipes, Section 13.8)

> Since 70 points are contained within a period one understands that 5000
> points approximatively correspond to 71 full periods. Note however, and
> this seems to be the main difficulty, that the time abscissas of the
> points are not equally spaced. Otherwise discrete Fourier transform
> should be convenient.

The Lomb periodogram is designed for this situation.

> Of course I could interpolate the function, rmod versus t, but the
> accuracy of the period obtained in that way is ridicoulusly low compared
> to the, say 50 figures, injected in the data.

Why?  Interpolation[data] generates an InterpolatingFunction object which
returns values with the same precision as those in data.

I've appended below a Notebook which gives an example of the
Lomb periodogram approach.

Cheers,
    Paul

____________________________________________________________________
Paul Abbott                                   Phone: +61-8-9380-2734
Department of Physics                           Fax: +61-8-9380-1014
The University of Western Australia
Nedlands WA  6907                     mailto:paul at physics.uwa.edu.au
AUSTRALIA                            http://physics.uwa.edu.au/~paul

            God IS a weakly left-handed dice player
____________________________________________________________________

Notebook[{
Cell["For the period:", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`\[Tau] = 1.2`20\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`1.2`20\)], "Output"]
}, Open  ]],

Cell["consider the following irregularly sampled (periodic) data:",
"Text"],

Cell[BoxData[
    \(TraditionalForm\`\(data =
        Function[{t}, {t, \(1 + 2 \( cos(t\/\[Tau])\)\)\/\(2 +
sin(t\/\[Tau])\
\)}] /@ SetPrecision[Sort[Table[Random[Real, {0, 50}], {500}]],
            20];\)\)], "Input"],

Cell[CellGroupData[{

Cell["Visualization", "Section"],

Cell["It is clear from the plot that the data is periodic.", "Text"],

Cell[BoxData[
    \(TraditionalForm\`\(ListPlot[data,
        PlotStyle -> AbsolutePointSize[1]];\)\)], "Input"],

Cell["Here we compare with the exact functional form.", "Text"],

Cell[BoxData[
    \(TraditionalForm\`\(dataplot =
        Plot[\(1 + 2 \( cos(t\/\[Tau])\)\)\/\(2 + sin(t\/\[Tau])\), {t, 0,
            50}, PlotStyle -> {AbsoluteThickness[1], Hue[0]},
          PlotPoints -> 30, Epilog -> First[%]];\)\)], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Spectral Energy", "Section"],

Cell[TextData[{
  "It is straightforward to compute the spectral energy \[LongDash]",
  " which is closely related to the ",
  "Lomb ",
  StyleBox["normalized periodogram",
    FontSlant->"Italic"],
  " (see ",
  StyleBox["Numerical Recipes",
    FontSlant->"Italic"],
  " \[Section]13.8 and ",
  StyleBox["The Mathematica Journal ",
    FontSlant->"Italic"],
  StyleBox["7",
    FontWeight->"Bold"],
  "(1):27-28) \[LongDash] using ",
  StyleBox["integration by pattern-matching",
    FontSlant->"Italic"],
  ":"
}], "Text",
  CellTags->"spectral energy"],

Cell[BoxData[
    \(TraditionalForm\`\[ScriptCapitalE](\[Omega]_) :=
      Fit[data, {sin(\[Omega]\ t), cos(\[Omega]\ t)},
          t] /. \[InvisibleSpace]\(cos(t_)\)\ \[Beta]_ + \[Alpha]_\ \(sin(
                t_)\) :> 1\/2\ \((\[Alpha]\^2 + \[Beta]\^2)\)\)], "Input"],

Cell["\<\
and then plot it as a function of \[Omega].\
\>", "Text"],

Cell[BoxData[
    \(TraditionalForm\`\(Plot[\[ScriptCapitalE](\[Omega]), {\[Omega], 0.1,
          2}, PlotRange -> All, AxesOrigin -> {0, 0}];\)\)], "Input"],

Cell["\<\
We can see the dominant frequency (and harmonics) and it is \
straightforward to determine the period:\
\>", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`1/\[Omega] /. \[InvisibleSpace]Last[
        FindMinimum[\(-\(\[ScriptCapitalE](\[Omega])\)\), {\[Omega], 0.7,
            0.75}]]\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`1.1979772479559643`\)], "Output"]
}, Open  ]],

Cell["\<\
However, even with such precise data, this method does not \
(cannot?) give an accurate value for the period. The fractional error \
is\
\>", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
    \(TraditionalForm\`\(% - \[Tau]\)\/\[Tau]\)], "Input"],

Cell[BoxData[
    \(TraditionalForm\`\(-0.0016856267033630672`\)\)], "Output"]
}, Open  ]],

Cell["which is wrong by over 0.16%.", "Text"]
}, Closed]]
}
]





  • Prev by Date: Re: searchforperiod
  • Next by Date: Re: Where's the Speed?
  • Previous by thread: Re: searchforperiod
  • Next by thread: Re: searchforperiod