6
$\begingroup$

I have irregularly sampled data that looks as in the picture.

Is there a function in Mathematica that allows me to find an approximate or average period or wavelength?

I know beforehand that the data will always look like this, with some random wiggling.

If there is no function that does this, how best to do it in Mathematica?

data - note the sampling rate increases with time

This is a link to the data. The data is very sparse at the beginning and grows denser towards the end. Use ListLinePlot rather than ListPlot to visualize:

https://www.wolframcloud.com/obj/egarcia/Published/M%20STACK%20DATA%202022%20may%2025.nb

data=Import["https://pastebin.com/raw/JvUgqjxT"]

EDIT It is hard to choose a best solution. Perhaps D. Lichtbau's solution using the irregular periodogram function in the function repository I like the most, but I also very much like the fitting solution, as it returns all the parameters needed to make a visual check of the hypothesized periodicity of the data. And as remarked by Roman in a comment, perhaps the best is to use these two in combination. For the sake of just choosing one, I will choose Daniel Lichtbau's solution, perhaps because it is the closest thing to having a function within the Wolfram Language that does what we want.

$\endgroup$
8
  • 4
    $\begingroup$ It would be a good idea to post the data. If too long, place it somewhere accessible and post a link. I for one would not try to tackle this from a picture/image alone. $\endgroup$ Commented May 24, 2022 at 22:46
  • 3
    $\begingroup$ The obvious approach would be to use the Fourier command. $\endgroup$
    – bill s
    Commented May 25, 2022 at 0:12
  • 2
    $\begingroup$ Reading this may be important, it's nit that simple to address the non-uniform sampling $\endgroup$
    – rhermans
    Commented May 25, 2022 at 8:43
  • 2
    $\begingroup$ @EGME I have edited your question to provide an easy way for other users to access your data. In this case an Import function with the address of a permanent repository of data (Pastebin). In the future, please do make an effort to facilitate the work for people trying to help you. $\endgroup$
    – rhermans
    Commented May 25, 2022 at 9:02
  • 2
    $\begingroup$ Also there are prior MSE posts involving FT-like computations on irregularly spaced data. The one here contains several salient features insofar as responses cover a few ways to go about this. $\endgroup$ Commented May 25, 2022 at 23:02

4 Answers 4

7
$\begingroup$

This is something one can tackle using an irregular periodogram, and that happens to be available in the Wolfram Function Repository.

Grab the data and separate into times and values.

data = Import["https://pastebin.com/raw/JvUgqjxT"];
{times, vals} = Transpose[data];

Plot a periodogram. From picture in the post, guessing we should go to frequency as high as 25.

pgram = Plot[
  ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, 0., 
   25}, PlotRange -> All]

enter image description here

So we see a peak in the vicinity of 14. We can home in on that.

{max, freq} = 
 FindMaximum[
  ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, 14}]

(* Out[23]= {0.0284429, {w -> 14.0784}} *)

the approximated period can be computed from this.

per = 2*Pi/w /. freq

(* Out[24]= 0.4463 *)

Another way to approximate the period, directly in this case, is to use the Lafler-Kinman (or string length) method.

pgramLK = 
 Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals, 
   Method -> "LaflerKinman"], {w, 0., 1}, PlotRange -> All]

enter image description here

Now we can minimize this to approximate the period. Here I use NMinimize to get around a problematic gradient computation. As there are two dips (indicating we might have to halve the main frequency) we'll compute both.

{min, per} = 
 NMinimize[
  ResourceFunction["IrregularPeriodogram"][w, times, vals, 
   Method -> "LaflerKinman"], {w, .4, .45}]

(* Out[34]= {5.17042, {w -> 0.447938}} *)

{min2, per2} = 
 NMinimize[
  ResourceFunction["IrregularPeriodogram"][w, times, vals, 
   Method -> "LaflerKinman"], {w, .85, .95}]

(* Out[43]= {3.62404, {w -> 0.888249}} *)

Check that they are approximately consistent.

w/2 /. per2

(* Out[45]= 0.444125 *)

I will remark that I trust more the default method, in the sense that it gives no indication of lower frequency peak near 7.07.

$\endgroup$
14
  • $\begingroup$ Thank you. This seems quite very useful, and perhaps the best solution yet. $\endgroup$
    – EGME
    Commented May 25, 2022 at 15:41
  • 1
    $\begingroup$ I was about to try to implement Lomb-Scargle from scratch and then I saw this. Very nice. $\endgroup$ Commented May 25, 2022 at 18:28
  • $\begingroup$ Hello Daniel, I cannot get my computer to produce the plot of the periodogram as in your first figure. Is there a special incantation that is needed for that to work? Sorry if this is a bad question. $\endgroup$
    – EGME
    Commented May 26, 2022 at 17:06
  • $\begingroup$ The question is fine. Did you follow all steps I showed? If so, what version are you using? It is possible that I have version dependent code in that function (though offhand I don't recall any such). $\endgroup$ Commented May 26, 2022 at 18:15
  • $\begingroup$ Yes, I followed all the steps. I am using version 13.0.1, the latest. It computes the periodogram all right, I get the value I need. But I cannot plot it. $\endgroup$
    – EGME
    Commented May 26, 2022 at 21:03
5
$\begingroup$

Maybe a nonlinear fit would work for you, but it needs a good starting point to converge:

f = NonlinearModelFit[data,
                      a + b Sin[c x + d],
                      {{a, 0}, {b, 0.01}, {c, 14}, {d, 2.5}},
                      x];
f["BestFitParameters"]
(*    {a -> -0.00119216, b -> 0.0105601, c -> 14.0715, d -> 2.47229}    *)

Plot[f[x], {x, Sequence @@ MinMax[data[[All, 1]]]}, Epilog -> {Red, Point[data]}]

enter image description here

$\endgroup$
4
  • $\begingroup$ Thank you, it looks like this might work. I will try it with other datasets when they become available. However, how do you obtain the "period" with your method? In theory, this should be almost periodical. $\endgroup$
    – EGME
    Commented May 25, 2022 at 8:05
  • $\begingroup$ Ok, right, sorry, that info is provided. I will accept this answer after 24 hours. I would like to see if there are any other suggestions. $\endgroup$
    – EGME
    Commented May 25, 2022 at 8:26
  • $\begingroup$ You can get the period with 2π/c /. f["BestFitParameters"]. $\endgroup$
    – Roman
    Commented May 25, 2022 at 9:08
  • $\begingroup$ Thanks, I realized it, but I was unable to edit the comment because more than 5 minutes had passed! $\endgroup$
    – EGME
    Commented May 25, 2022 at 9:15
4
$\begingroup$

Try using Periodogram. See for more HERE. In short on an example, let's take temperature data - about 2 decades, sampling monthly:

temp=WeatherData["Chicago","Temperature",{{2000,1,1},{2021,6,31},"Month"}];

You got something pretty periodic, but not exactly periodic:

DateListPlot[temp,FrameLabel->Automatic,Mesh->All,MeshStyle->Red]

enter image description here

Now Periodogram picks frequency related to 1-year period easily with a peak at 1:

Periodogram[QuantityMagnitude[Values@temp],
SampleRate->12,
PlotRange->All,Frame -> True,AspectRatio -> 1/3,
GridLines -> {{1}, None}, 
GridLinesStyle -> Directive[Red,Thick, Dashed]]

enter image description here

Ignore first peak at zero, it is probably related to the fact that the data have none-zero average. To understand why this works option SampleRate->12 is important. Your data basically are something like transformation of a basic function similar to

Sin[2 Pi / (1 year)]

so your frequency is 1 and you sample with a month meaning 1/12-th of a year.

$\endgroup$
1
  • $\begingroup$ I don't think your idea will work because the sampling rate is not fixed. It increases (exponentially?) with time. I changed the picture and added a link to a notebook with the data $\endgroup$
    – EGME
    Commented May 25, 2022 at 6:54
2
$\begingroup$

To get an equal distant sampling you can use Interpolation and resample (assuming your data is stored in "data" and I eliminate the x values because they are not important for frequencies for equidistant data points):

dat = Table[int[x], {x, 10.05, 13.8, 0.01}];
ListPlot[dat]

enter image description here

If we now get the periodogram we see that there is no dominant frequency.

Periodogram[dat, PlotRange -> All] 

enter image description here

$\endgroup$
10
  • $\begingroup$ Could you elaborate on how to use Interpolation to resample? Also, which would be the lowest frequency in the Periodogram? I am not familiar with this function. I scanned the documentation, but I will have to study what its output means. What I want is the lowest frequency (not a dominant one). I think this is a neat solution if I can get from it the information that I want. For example, in one of the solutions above you get a parameter c=14.075 from which you can get the period. If I can get that information with your method, it would be almost perfect. $\endgroup$
    – EGME
    Commented May 25, 2022 at 9:23
  • $\begingroup$ I also think that if you have a lot of data, your solution would be faster to compute. Fitting can take a long time. $\endgroup$
    – EGME
    Commented May 25, 2022 at 9:23
  • $\begingroup$ Right, I see how to use Interpolation now. So if I can get the period from the periodogram, that would be quite a very good solution. In the x-coordinate scale of your plot above would be fine, I can translate that into the original scale, I think $\endgroup$
    – EGME
    Commented May 25, 2022 at 9:34
  • $\begingroup$ The x coordinates goes from 0 to 0.5. You need at least 2=1/0.5 points in a period to sample a frequency (Nyquist theorem). The lowest frequency is approx. 0.02, that gives therefore a period of approx. 1/0.02=50 data points. $\endgroup$ Commented May 25, 2022 at 9:41
  • $\begingroup$ Note also, zero frequency is the DC part of the signal (the mean). $\endgroup$ Commented May 25, 2022 at 9:43

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.