Wolfram Research

Function Repository Resource:

IrregularPeriodogram

Source Notebook

Compute a periodogram for data from unevenly spaced intervals

Contributed by: Daniel Lichtblau

ResourceFunction["IrregularPeriodogram"][w,times,vals]

computes the irregular periodogram at frequency w for the given set of times and vals.

ResourceFunction["IrregularPeriodogram"][w,mat]

computes the irregular periodogram for the set of time and value pairs in the two-column matrix mat.

Details and Options

IrregularPeriodogram accepts a Method option.
Allowed Method values are "Deeming", "LombScargle", "LaflerKinman", and "StringLength", where the last two are equivalent.
The Deeming periodogram for times {t1,,t1} and corresponding values {v1,,v1} is given by 𝒫(ω)=(∑j=1nvj-ⅈ ω tj❘)2.
Unlike the classical case, the Lomb-Scargle periodogram here does not normalize by the variance of the data values.
Strictly speaking, the "LaflerKinman" setting will not give a periodogram computation. Instead of computing a value based on a given frequency, it computes directly based on a putative period p. Small values indicate a short total length when folding by p along the time axis (that is, sorting the times modulo p and approximating the total length between consecutive time-folded points).

Examples

Basic Examples

Take random times and corresponding values for the function f(t)=sin(t)+cos(3t)/2:

In[1]:=
SeedRandom[11222333344444]
data = Sort[
   Map[{#, Sin[#] + Cos[3*#]/2} &, RandomReal[{0, 20}, 500]]];

Subtract the mean of the data values and plot the resulting time series:

In[2]:=
times = data[[All, 1]];
vals = data[[All, 2]] - Mean[data[[All, 2]]];
ListPlot[Transpose[{times, vals}]]
Out[4]=

Plot the periodogram computed from this unevenly spaced set of measurements:

In[5]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, .1,
   5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[5]=

Compute the larger frequency (the location of the larger spike):

In[6]:=
freq = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals], .7 <= w <= 1.3}, w]
Out[6]=

Compute the period from this frequency:

In[7]:=
per = 2 Pi/freq
Out[7]=

Plot the time series “folded” by this period (times are computed modulo the computed period) to see the approximately sinusoidal variation in the light curve:

In[8]:=
foldTimes[times_, period_] := period*FractionalPart[times/period]
In[9]:=
ListPlot[Transpose[{foldTimes[times, per], vals}]]
Out[9]=

The string-length approach also shows this periodicity:

In[10]:=
Plot[ResourceFunction["IrregularPeriodogram"][p, times, vals, Method -> "LaflerKinman"], {p, 1, 8}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"p", "\!\(\*SubscriptBox[\(P\), \(LK\)]\)(p)"}, ImageSize -> 400]
Out[10]=

One can get a fairly good estimate of the frequencies present using far fewer samples, so long as they are irregularly spaced in time. We take 20 random points from the original set above:

In[11]:=
SeedRandom[11222333344444]
data2 = Sort[RandomSample[data, 20]];
times2 = data2[[All, 1]];
vals2 = data2[[All, 2]] - Mean[data2[[All, 2]]];
ListPlot[Transpose[{times2, vals2}]]
Out[12]=

Plot the periodogram from this subset:

In[13]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times2, vals2], {w, .1, 5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[13]=

Estimate the period for the larger frequency:

In[14]:=
freq = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times2, vals2], .7 <= w <= 1.3}, w]
Out[14]=

Scope

We use a light curve time series from the Wolfram Demonstration “CepheidVariableStarLightCurve”:

Subtract the mean of the data values and plot the resulting time series:

In[15]:=
times = data[[All, 1]];
vals = data[[All, 2]] - Mean[data[[All, 2]]];
ListPlot[Transpose[{times, vals}]]
Out[17]=

Plot the periodogram computed from this unevenly spaced set of measurements:

In[18]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, .1,
   5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/days)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[18]=

Compute the main frequency as the location of the larger spike:

In[19]:=
freq = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals], 2. <= w <= 3}, w]
Out[19]=

Compute the period from this frequency:

In[20]:=
per = 2 Pi/freq
Out[20]=

Plot the time series “folded” by this period (times are computed modulo the computed period) to see the approximately sinusoidal variation in the light curve:

In[21]:=
ListPlot[Transpose[{foldTimes[times, per], vals}]]
Out[21]=

Plot the Lomb-Scargle periodogram for the same data:

In[22]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals, Method -> "LombScargle"], {w, .1, 5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/days)", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[22]=

The computed main frequency and period agree with the values from the Deeming periodogram:

In[23]:=
freqLS = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals, Method -> "LombScargle"], 2.5 <= w <= 3}, w];
perLS = 2 Pi/freqLS;
{freqLS, perLS}
Out[24]=

Plot the string-length values for the same data:

In[25]:=
Plot[ResourceFunction["IrregularPeriodogram"][p, times, vals, Method -> "LaflerKinman"], {p, .1, 5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"p (days)", "\!\(\*SubscriptBox[\(L\), \(LK\)]\)(p)"}, ImageSize -> 400]
Out[25]=

Compute the best estimate for the period based on the minimum value attained by the Lafler-Kinman computation for this data:

In[26]:=
perLK = NArgMin[{ResourceFunction["IrregularPeriodogram"][p, times, vals, Method -> "LaflerKinman"], 2. <= p <= 2.5}, p]
Out[26]=

A periodogram can often detect multiple frequencies provided they are well separated, even if they are not commensurate (that is, there is no actual “period” for the data). We extract values from such a function with frequencies at :

In[27]:=
SeedRandom[11222333344444]
data = Sort[
   Map[{#, Sin[#] - 2*Sin[E*#] Cos[Sqrt[2]*#] + Cos[Pi*#]/2} &, RandomReal[{0, 20}, 800]]];
times = data[[All, 1]];
vals = data[[All, 2]] - Mean[data[[All, 2]]];
ListPlot[Transpose[{times, vals}]]
Out[28]=

Plot the periodogram for this data:

In[29]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, .1,
   5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[29]=

We get one estimated frequency around midway between 1 and :

In[30]:=
NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals], .7 <= w <= 1.4}, w]
Out[30]=

Another is approximately π:

In[31]:=
NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals], 3 <= w <= 3.3}, w]
Out[31]=

A third is close to :

In[32]:=
NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals], 3.9 <= w <= 4.6}, w]
Out[32]=

We can get rough estimates with far fewer sampled values:

In[33]:=
SeedRandom[11222333344444]
data = Sort[
   Map[{#, Sin[#] - 2*Sin[E*#] Cos[Sqrt[2]*#] + Cos[Pi*#]/2} &, RandomReal[{0, 20}, 20]]];
times = data[[All, 1]];
vals = data[[All, 2]] - Mean[data[[All, 2]]];
ListPlot[Transpose[{times, vals}]]
Out[34]=

Peaks of the periodogram give estimates of the frequencies:

In[35]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, .1,
   5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[35]=

Properties & Relations

When the data is equally spaced on the horizontal axis, the irregular periodogram is essentially the square of the Fourier transform; we modify a previous example to demonstrate this:

In[36]:=
SeedRandom[11222333344444]
times = Range[2, 20, .1];
data = Table[Sin[j] + Cos[3*j]/2 + .01*RandomReal[], {j, times}];

Subtract the mean of the data values and plot the resulting time series:

In[37]:=
vals = data - Mean[data];
ListPlot[vals]
Out[38]=

Plot the periodogram computed from this (evenly) spaced set of measurements:

In[39]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals], {w, .1,
   5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[39]=

The square of the Fourier transform gives a similar result once we properly scale as frequencies on the horizontal axis:

In[40]:=
ftsqr = Abs[Fourier[vals]]^2;
ListPlot[ftsqr[[1 ;; 20]], DataRange -> {0, 2 Pi}, PlotRange -> All]
Out[41]=

In the Wolfram Language this can also be computed using PeriodogramArray:

In[42]:=
ListPlot[PeriodogramArray[vals][[1 ;; 20]], DataRange -> {0, 2 Pi}, PlotRange -> All]
Out[42]=

Yet another way to obtain these absolute values is with the Fourier transform of the set of values convolved with itself:

In[43]:=
ftconvolved = Fourier[ListConvolve[vals, vals, {1, 1}, vals]]/Sqrt[Length[vals]];
ListPlot[Abs@ftconvolved[[1 ;; 20]], DataRange -> {0, 2 Pi}, PlotRange -> All]
Out[44]=

The Lafler-Kinman string-lengths-based measure can be made to more closely resemble an ordinary periodogram using reciprocals. We use a basic example:

In[45]:=
SeedRandom[11222333344444]
data = Sort[Map[{#, Sin[#] + Cos[3*#]/2} &, RandomReal[{0, 20}, 500]]];
times = data[[All, 1]];
vals = data[[All, 2]] - Mean[data[[All, 2]]];
Plot[irregPG[w, times, vals], {w, .1, 5}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}]
Out[46]=

Now again show the lafler-Kinman plot using this data:

In[47]:=
Plot[irregPG[p, times, vals, Method -> "LaflerKinman"], {p, 1, 8}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"p", "\!\(\*SubscriptBox[\(P\), \(LK\)]\)(p)"}]
Out[47]=

We go from time to frequency using the formula f=2π/t, and likewise go from small to large on the vertical axis by taking reciprocals:

In[48]:=
FrequencyAmplitude[{time_, len_}] := {2*Pi/time, 1/len}

Generate and plot a list of {foldtime,magnitude} pairs:

In[49]:=
pairs = Table[{p, irregPG[p, times, vals, Method -> "LaflerKinman"]}, {p, 1., 15, 1/400}];
ListPlot[pairs]
Out[50]=

Now do the transformation and replot:

In[51]:=
ListPlot[Map[FrequencyAmplitude, pairs], PlotRange -> All]
Out[51]=

We see that the main spike coincides with that of the Deeming periodogram. We also see a small bump in the Deeming periodogram near frequency 3. For the Lafler-Kinman version this corresponds to a fold time near 2.1:

In[52]:=
NSolve[2. Pi/t == 3]
Out[52]=

Compute a lot of values in the vicinity of this floding time:

In[53]:=
pairsNearby = Table[{p, irregPG[p, times, vals, Method -> "LaflerKinman"]}, {p, 1.85, 2.3, 1/10000}];
ListPlot[Map[FrequencyAmplitude, pairsNearby], PlotRange -> All]
Out[54]=

Applications

Sunspot activity cycles

Sunspot activity follows an approximately 11 year cycle. Download a century of monthly averaged values from the Wolfram Data Repository:

In[55]:=
monthsunspots = ResourceFunction["SpaceWeatherData"][{DateObject[{1918}], DateObject[{2018}]}, "MonthlySunspotNumber"];

Show the plot for this time series:

In[56]:=
DateListPlot[monthsunspots]
Out[56]=

Obtain the times and mean-centered values:

In[57]:=
tvpairs = monthsunspots["Path"];
tlist = tvpairs[[All, 1]];
vals = tvpairs[[All, 2]] - Mean[tvpairs[[All, 2]]];
times = N[(tlist - First[tlist])/(tlist[[2]] - First[tlist])];

Plot the Deeming periodogram for this time series:

In[58]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals, Method -> "Deeming"], {w, .001, .1}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/months)", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[58]=

Compute the maximum frequency and the corresponding estimated period in months:

In[59]:=
freqD = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals, Method -> "Deeming"], .045 <= w <= .06}, w];
perD = 2*Pi/freqD;
{freqD, perD}
Out[60]=

Plot the Lomb-Scargle periodogram for this time series:

In[61]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times, vals, Method -> "LombScargle"], {w, .001, .1}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/months)", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[61]=

Compute the maximum frequency and the corresponding estimated period in months:

In[62]:=
freqLS = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times, vals, Method -> "LombScargle"], .04 <= w <= .06}, w];
perLS = 2*Pi/freqLS;
{freqLS, perLS}
Out[63]=

Plot a list of string-lengths vs. folding period for this time series:

In[64]:=
ListPlot[lkvals = Table[ResourceFunction["IrregularPeriodogram"][p, times, vals, Method -> "LaflerKinman"], {p, 1, 300, .1}]]
Out[64]=

Compute the estimated period from the low point position:

In[65]:=
perLK = 1 + (FirstPosition[lkvals, Min[lkvals]][[1]] - 1)*.1
Out[65]=

All three estimates are in the ballpark of 10.5 years. We plot the time series epoch-folded by the second and third period estimates, as well as by the conventional 11 years:

In[66]:=
GraphicsRow[{ListPlot[Transpose[{foldTimes[times, perLS], vals}]], ListPlot[Transpose[{foldTimes[times, perLK], vals}]], ListPlot[Transpose[{foldTimes[times, 12*11], vals}]]}]
Out[66]=

Another cepheid variable star analysis

Download and process times and values of the light curve for the variable star known as cep2308:

In[67]:=
cep2308Idata = Import["http://ogledb.astrouw.edu.pl/~ogle/CVS/data/I/08/OGLE-LMC-\
CEP-2308.dat", "Table"][[All, 1 ;; 2]];
{times2308, ovalues2308} = Transpose[cep2308Idata];
times2308 = times2308 - times2308[[1]];
mean2308 = Mean[ovalues2308];
values2308 = ovalues2308 - mean2308;

Show a plot of measured amplitudes over time:

ListPlot[Transpose[{times2308, values2308}],AxesLabel{"time (days)","amplitude fluctuation"},ImageSize400]
Out[68]=

Plot the Deeming periodogram and compute the corresponding estimated period:

In[69]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times2308, values2308, Method -> "Deeming"], {w, .001, 3}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/days)", "\!\(\*SubscriptBox[\(P\), \(D\)]\)(\[Omega])"}, ImageSize -> 400]
Out[69]=
In[70]:=
freqD = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times2308, values2308], .8 <= w <= 1.3}, w];
perD = 2 Pi/freqD;
{freqD, perD}
Out[70]=

Plot the Lomb-Scargle periodogram and compute the corresponding estimated period:

In[71]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times2308, values2308, Method -> "LombScargle"], {w, .001, 3}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/days)", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[71]=
In[72]:=
freqLS = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, times2308, values2308, Method -> "LombScargle"], .8 <= w <= 1.3},
    w];
perLS = 2 Pi/freqLS;
{freqLS, perLS}
Out[72]=

Plot the Lafler-Kinman folded epoch total string lengths and compute the corresponding estimated period:

In[73]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, times2308, values2308, Method -> "LaflerKinman"], {w, .001, 20}, PlotPoints -> 300, PlotRange -> All, AxesLabel -> {"\[Omega] (radians/months)", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[73]=
In[74]:=
perLK = NArgMin[{ResourceFunction["IrregularPeriodogram"][p, times2308, values2308, Method -> "LaflerKinman"], 5.4 <= p <= 6}, p]
Out[74]=

The period estimates all agree to several places. We plot the time series folded by epoch using the second one:

In[75]:=
ListPlot[Transpose[{foldTimes[times2308, perLS] , values2308}], AxesLabel -> {"time (days)", "amplitude fluctuation"}, ImageSize -> 400]
Out[75]=

Searching for periodicity in a chromosome

Coding sections of genes are well known to exhibit a periodicity of 3 base pairs (bp). Often there are other periodicities. We show a method for finding them, using chromosome 12 from Saccharomyces cerevisiae (baker’s yeast) as an example:

In[76]:=
seq = Import[
   "ftp://ftp.ensembl.org/pub/release-81/fasta/saccharomyces_\
cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna_rm.chromosome.XII.\
fa.gz"];
chromo = seq[[1]];

Literature on the subject states that it has a periodicity in the 10-11 range, when looking at tetramers of the form A*T*, where the star here denotes zero or more occurrances of the preceding nucleotide (some specific values reported are 10, 10.2, 10.4, and 10.5). Split the chromosome into segments of 1000 bp, and locate positions of the 4-mers “AAAA”, ”AAAT”, ”AATT”, ”ATTT” and ”TTTT”:

In[77]:=
chromos = StringSplit[chromo, "N"] /. "" :> Nothing;
sublengths = 1000;
subchromos = Flatten[Map[StringPartition[#, sublengths] &, chromos]];
at4 = "AAAA" | "AAAT" | "AATT" | "ATTT" | "TTTT";
posnlists = Map[StringPosition[#, at4][[All, 1]] &, subchromos];

While there are various ways to turn a gene sequence into a numerical sequence, we instead employ the irregular periodogram on common gaps between the positions we found:

In[78]:=
pdiffs = Flatten[
   Table[Outer[Subtract, posnsj, posnsj], {posnsj, posnlists}]];
pdiffs = DeleteCases[pdiffs, aa_ /; aa <= 0];
tallies = Tally[pdiffs];
talliesS = Sort[tallies];

We discard smaller counts if there are more than 500. We further decimate by removing all gap counts that are smaller than at least one of their 19 subsequent neighbors:

In[79]:=
grouped = Partition[talliesS, 20, 1];
upper = 500;
groupedS = DeleteDuplicates[
   Map[First[SortBy[#, -#[[2]] &]] &, grouped[[2 ;; UpTo[upper]]]]];

Plot number of gaps as a function of gap size:

In[80]:=
ListPlot[groupedS, AxesLabel -> {"gap", "bin size"}, ImageSize -> 400]
Out[80]=

Treat gap lengths as the “time” variable and gap counts as the dependent variable:

In[81]:=
{timesC, valsC} = Transpose[groupedS];
valsC = valsC - Mean[N@valsC];

We already know of a periodicity of 3 and seek larger ones, so we plot up to a frequency range that will only give periods of at least 4:

In[82]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, timesC, valsC, Method -> "LombScargle"], {w, .3, 2*Pi/4}, PlotPoints -> 500, PlotRange -> All, AxesLabel -> {"\[Omega]", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[82]=

Isolate the large frequency near 0.6 to find a plausible periodicity value:

In[83]:=
freqLS = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, timesC, valsC, Method -> "LombScargle"], .5 <= w <= .7}, w];
perLS = 2*Pi/freqLS;
{freqLS, perLS}
Out[84]=

Homing in near 0.6 in the plot will show that there is a double peak, with the slightly larger one giving rise to a periodicity less than 10. This may be due to an interaction between the periodicity of 3 and the larger one. We can attempt to curtail this effect by repeating the periodogram computations on data that has had gaps that are multiples of 3 removed:

In[85]:=
{timesC, valsC} = Transpose[Select[groupedS, (! IntegerQ[#[[1]]/3]) &]];
valsC = valsC - Mean[N@valsC];
Plot[ResourceFunction["IrregularPeriodogram"][w, timesC, valsC, Method -> "LombScargle"], {w, .3, 2*Pi/4}, PlotPoints -> 500, PlotRange -> All, AxesLabel -> {"\[Omega]", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[87]=

Home in near 0.6:

In[88]:=
Plot[ResourceFunction["IrregularPeriodogram"][w, timesC, valsC, Method -> "LombScargle"], {w, .4, .8}, PlotPoints -> 500, PlotRange -> All, AxesLabel -> {"\[Omega]", "\!\(\*SubscriptBox[\(P\), \(LS\)]\)(\[Omega])"}, ImageSize -> 400]
Out[88]=

Recompute the period estimate:

In[89]:=
freqLS = NArgMax[{ResourceFunction["IrregularPeriodogram"][w, timesC, valsC, Method -> "LombScargle"], .4 <= w <= .8}, w];
perLS = 2*Pi/freqLS;
{freqLS, perLS}
Out[91]=

This is in line with values reported in the literature. A string-length plot and optimization gives a value in the same vicinity:

In[92]:=
Plot[ResourceFunction["IrregularPeriodogram"][p, timesC, valsC, Method -> "LaflerKinman"], {p, 8, 11}, PlotRange -> All]
Out[92]=
In[93]:=
NArgMin[{ResourceFunction["IrregularPeriodogram"][p, timesC, valsC, Method -> "LaflerKinman"], 10 <= p <= 11}, p]
Out[93]=

Different settings will give slightly different estimates, for example, with sublengths set to 5000 and upper set to 2000 the estimate is around 10.28. So the method here is viable but only to low precision.

Resource History

License Information