Allan Deviation Algorithm

Question

I am trying to write my own Allan Deviation calculator in Mathematica using the definition, $$\sigma(\tau)^{2}=\frac{1}{2M-2}\sum_{i=1}^{N-1}(\bar{y}_{i+1} - \bar{y}_{i})^{2}\text{.}$$ Where the $\bar{y}$ are mean partitioned data for some multiple of measurement time $\tau$, and $M$ is the length of the number of $(\bar{y}_{i+1} - \bar{y}_{i})^{2}$ associated with each $\tau$. I made a toy model of partitioning the data that I am happy with, but I am struggling to extending it to the next step, which is simply summing the $(\bar{y}_{i+1} - \bar{y}_{i})^{2}$ for each $\tau$. The toy model is:

NumberPoints = 64;
IndexData = Range [NumberPoints];

Do[
  Print["________________τ=", j];
  Do[
   Print[
    "(MEAN", Partition[IndexData, j][[2 i]], "-", "MEAN", 
    Partition[IndexData, j][[2 i - 1]], 
    "\!\(\*SuperscriptBox[\()\), \(2\)]\)"
    ],
   {i, 1, (NumberPoints/(j*2))}
   ],
  {j, 1, NumberPoints/4}
  ];

If I try, for the next step:

Do[
  Clear[BlockMeanDifferenceList];
  BlockMeanDifferenceList = {};
  Do[
   BlockMeanDifferenceList = (
     Mean[Partition[RandomData, τ][[2 i]]] - 
     Mean[Partition[RandomData, τ][[2 i - 1]]]
   )^2;
   AppendTo[σ, (1/(2 Length[BlockMeanDifference] - 2)*Total[BlockMeanDifference])^(1/2)];,
   {i, 1, (NumberOfPoints/(τ*2))}
   ];,
  {τ, 1, NumberOfPoints/4}
  ];

My sigma list is returned for zero in all elements...

Any thoughts?

Answer 1

Below an attempt to implement your formula as I understand it (the prefactor is missing, as I'm not quite sure what $M$ should be) - please do correct me if I misunderstood something.

avar[data_, τ_] := Total[
  (
    Subtract @@@ (*get the difference of the pairs*)
     Partition[ (*create overlapping subsequences of length 2*)
      Mean /@ Partition[data, τ], (*partition data into block and calculate mean*)
      2, 1
      ]
    )^2
  ]

With this, avar[Array[Subscript[a, #] &, 16], 4] yields:

If you don't want overlapping differences, simply remove the ,1 in the outer partition (this controls the offset of the individual sublists)

Update

Since you tagged this procedural-programming, here is an approach more closely based on yours (it's impossible to tell you what's wrong with your attempt specifically, as there are still a few obvious typos - see my comment)

avarProc[data_, τ_] := Module[
  {
   partitions = Partition[data, τ],
   σ = {},
   diff
   },
  Do[
   diff = (Mean@partitions[[i + 1]] - Mean@partitions[[i]])^2;
   AppendTo[σ, diff],
   {i, 1, Length@data/τ - 1}
   ];
  Total@σ
  ]

This produces the same output as the above. Again, this uses overlapping sublists. To change this, simply take your indexing and loop bounds.