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]]
}
]