       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[
}, Open  ]],

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

Cell[BoxData[
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[
PlotStyle -> AbsolutePointSize];\)\)], "Input"],

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

Cell[BoxData[
Plot[\(1 + 2 \( cos(t\/\[Tau])\)\)\/\(2 + sin(t\/\[Tau])\), {t, 0,
50}, PlotStyle -> {AbsoluteThickness, Hue},
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[
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[
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[
FindMinimum[\(-\(\[ScriptCapitalE](\[Omega])\)\), {\[Omega], 0.7,
0.75}]]\)], "Input"],

Cell[BoxData[
}, 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[
}, 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