A Modest Proposal

* Texted Ray over the weekend about the possibility of designing a scintillator-based Geiger counter substitute that used the iPhone as a user interface. We talked about it on the phone for a little while this morning.

* Ray asked me to write up a paragraph or so of text for the proposal based on the recent results on the GPS timing system. Here's what I've got so far, about a page. Asked Ray for some feedback on what I should cut.

We recently completed a series of tests of the DeLorme GPS2058 module which we are using in our demonstration system for purposes of assigning absolute time stamps to individual shower detection events. We characterized the time stability of the module by comparing its output against a high-precision 10.0000 MHz OCXO (oven-controlled crystal oscillator) having 10 ppb (parts per billion) frequency stability, used as a frequency reference, and plotted the relative Allan deviation (a measure of frequency instability) between these two timing sources (below left). After correcting for the expected quantization error (due to the 10 Msps sampling rate) of ±50 ns, we found that the fine-grained Allan deviation (between adjacent 1-second frequency averaging windows) was 1.67 × 10⁻⁷ (in dimensionless units), which, at the 10 MHz OCXO frequency, translates into ±1.67 Hz of relative frequency instability, or in other words, ±167 ns RMS deviation from the mean in the time of arrival of individual PPS (pulse per second) timing pulses received from the GPS module. This then implies a comparable lower bound to the uncertainty of absolute time measurements based directly on the GPS PPS edges. We could improve the accuracy somewhat by averaging the GPS time readings over longer intervals, but (given the 10 MHz OCXO frequency) not to better than a ±25 ns lower bound that would be set by quantization error even if we used dual-edge triggered flip-flops in our input capture circuit to achieve 20 Msps. But without using differential GPS or military-grade access codes for improved accuracy, the GPS network does not offer users any better time resolution than around that level in any case.

If we assume, for the time being, an uncertainty of ∆t = ±160 ns for the absolute time measurements at a given detector site, then this limits the angular resolution for determinations of the bearing of a large-scale (atmospheric or interstellar) cosmic-ray shower detected at multiple sites. For two detectors separated by a straight baseline of length b < 7,926 km (diameter of Earth), and for an interstellar shower coming from a heading on the celestial sphere at an angle θ away from the baseline axis, the uncertainty in the measurement of θ is given by |∆θ| = c∆t/b (sin θ)⁻¹, in the approximation where
c∆t << b, and where the time reference at the other detector is assumed to be perfectly accurate. The same formula also applies to atmospheric showers for baselines that are less than that of a typical shower diameter at ground level (a few km). In the figure at above right, we plot this quantity for a variety of baseline lengths and source elevations (relative to the baseline axis). So for example, for two sites that are spaced 1,000 km apart, the angle θ for an interstellar shower coming from a direction nearly orthogonal to the baseline axis could be determined to within about 10 seconds of arc once the coincidence between the two detection events (occurring ~1 ms apart, or less) has been identified.

OK, now done a shorter (1 paragraph) description:

Our demonstration system for timestamping cosmic-ray shower events utilizes a DeLorme GPS2058 module for (once per second) absolute time synchronization, together with a Connor-Winfield 10.0000 MHz OCXO (oven-controlled crystal oscillator) with 10 ppb frequency stability for more fine-grained, higher-precision time measurements. We recently completed a series of tests on these components for purposes of verifying the absolute accuracy that we can expect to obtain for comparison of event times between detector sites; our empirical results and analysis will be described in detail in an IEEE journal article (currently in preparation). Our measurements of relative Allan deviation (frequency instability) are shown in the figure at left below. From this data, we inferred that there is an estimated ±167 ns RMS jitter in the time of arrival of individual timing pulses from the GPS unit, which is likely traceable to the low-offset phase noise of the TCXO (temperature-compensated crystal oscillator) which is used for timing purposes inside the GPS module. Given the implied time uncertainty of |∆t| ≈ 160 ns, we can calculate the uncertainty |∆θ| = c|∆t|/b (sin θ)⁻¹ in the inferred angle θ between a shower's heading and the baseline axis between two detector sites that are a distance b apart (below, right). So for example, given two detector sites spaced 1,000 km apart, the off-axis angle θ for an (interstellar) shower coming from a direction nearly orthogonal to the baseline axis could be determined to within about 10 seconds of arc, once the coincidence between the two detection events has been identified.

Also digging up some other graphs to post on project page... OK, all of the recent graphs (and some explanation discussion) are now posted at http://www.neutralino.org/projects/cosmici#TOC-Recent-Results.

Quantization Error Redux

Ray says that Sachin & Vaibhav are out of pocket for a little while. We'll wait and work on the board ourselves when we have time. In the meantime, I am reading a book on USB at home in the evenings, which will help us out if we need to debug or make changes in the communication.

My plan for today is to work on the new analysis of how the quantization error affects the Allan deviation by developing a program to do the calculations numerically.

I decided I wanted to first try just writing the new analysis in C/C++ (for faster development since I am most familiar with this language). I liked the NetBeans IDE that Gordon introduced me to, so I will try using it. I already reinstalled it with C/C++ support added a few weeks ago. Today, registered my installation with Oracle. Realized it needed Cygwin for the underlying gcc compiler. There is already a Cygwin environment in the Nios II IDE, but I decided instead of trying to use that one with NetBeans also, it would be cleaner to just install another one for general use. Did that, selecting gcc/make/standard libraries/gdb to install. Also included USB development support in case we need to do any raw USB coding in C. Created a new project "AnalyzeQuantErr"; now ready to start programming...

OK, I got to the core of my calculation, which is the following integral (expressed in Mathematica syntax):

Integrate[PDF[NormalDistribution[mu, sig], Phi2 - Phi1]*
PDF[NormalDistribution[mu, sig], Phi3 - Phi2], {Phi1, -(1/2), 1/2},
{Phi2, x2 - 1/2, x2 + 1/2}, {Phi3, x3 - 1/2, x3 + 1/2}]

Unfortunately, Mathematica is unable to integrate the result over Phi2, although it does manage to do the other two integrals. I tried all possible orders of integration, and did not achieve any better result. So, it looks as if I will have to numerically evaluate the following integral over Phi2:

Integrate[(1/4)*(Erf[(-(1/2) + mu - Phi2)/(Sqrt[2]*sig)] -
Erf[(1/2 + mu - Phi2)/(Sqrt[2]*sig)])*
(Erf[(-1 + 2*mu + 2*Phi2 - 2*x3)/(2*Sqrt[2]*sig)] -
Erf[(1 + 2*mu + 2*Phi2 - 2*x3)/(2*Sqrt[2]*sig)]),
{Phi2, -(1/2) + x2, 1/2 + x2}]

Fortunately, the error function erf() is included in the C math library already, so I don't have to write code to evaluate that as well.

OK, the program is done now, and I ran it and it produced sensible-looking results:

Here's the same graph on a log-log scale, which emphasizes the low end of the scale more:

As we can see, when the inherent frequency variation of the oscillator over the 1-second window size approaches 0, the measured Allan deviation due to quantization error approaches 1/2 cycle. This probably makes some kind of sense, since quantized phase measurements can be off by at most (+/-) 1/2 cycle, and 1/2 is also the RMS error of the quantized deviation that I obtained before by integrating the Dirac delta function under the assumption that the quantization errors of the 3 points x1,x2,x3 are independent & uniformly distributed in the range [-1/2,+1/2] cycle.

And, when the inherent frequency variation of the oscillator (for frequencies averaged over 1-second intervarls) increases well above +/- 1 cycle per second (standard deviation), the Allan deviation asymptotically approaches the actual standard deviation, which is correct since the quantization error then becomes relatively insignificant, and Allan deviation is designed to give the actual deviation in the case of a white noise model of the oscillator frequency fluctuations.

So, this seems like a good, solid theoretical result.

Meanwhile, let's apply this to our data set. We measured a 1-s Allan deviation of 1.7407 cycles. Going backwards through the above curve, we find that the actual standard deviation of frequency is 1.6674. Lo and behold, this is exactly the same as we had computed before by just using a simple quadrature combination of quantization error and inherent frequency variation, using the standard deviation of quantization error of 1/2 cycle that we had already obtained by our simpler analysis (integrating Dirac delta function)! This suggests that the above curve, arrived at by a complicated numerical integration procedure, could have been reproduced using the simple quadrature rule. Indeed, the curve obtained that way is indistinguishable to within 1/1000 cycle, which can be easily accounted for by the imprecision inherent to numerical integration.

This is all a good sanity check of the analysis - the two methods I used (simple quadrature rule & more detailed model) obtain the exact same result. So, I am now confident in this method (either/both method) for accounting for the quantization errors. Of course, it all assumes Gaussian white noise.

Deep Thought

Today (Wed., 3/30/11) I am taking notes and trying to figure out exactly how to *correctly* derive the expected Allan deviation given the quantization error and an oscillator model. I've decided I want to use a Gaussian white noise model for the inherent frequency fluctuations in the oscillator. (Other behavior should show up in the Allan deviations for larger window sizes.)

Key parameters are:

1. Standard deviation of oscillator frequency (when averaged over 1-s window sizes) in cycles per nominal second. (We can plot the result as a function of this, and use that to do a Bayesan back-calculation of the probabilities of different standard deviations given the raw deviation data. Actually, this dist would be sharply peaked, given the large amount of data available - so we can just do a direct reverse lookup of the standard deviation given the Allan deviation.)

2. Fractional part only of mean oscillator frequency, in [0,1] say, because the integer part is irrelevant to the quantization error.

And computational variables include:

3. The possible values of the quantized frequency deviation samples d1=x1-2x2+x3 for adjacent 1-second windows. These are small integers; a range [-10,+10] would be more than adequate. For this we scan through the range of possible values.

4. Given d1, always let x1=0, and then scan possible values of x2. These will always be integers in a range close to the integer part of mean oscillator frequency. For example, 10,000,012 +/- 10. Given d1 and x2, compute x3=d1+2*x2.

5. Phase error phi1 of first (x1) measurement in oscillator cycles, in the range [-1/2, 1/2]. This is uniformly distributed so we can just integrate over the range. The actual initial phase is Phi1=x1-phi1.

6. Actual frequencies F1=Phi2-Phi1 and F2=Phi3-Phi2 can then be integrated over the range of values that will yield the target values of x2 and x3. And, at each value, the probability density values can be looked up from the normal dist. Overall answer of the integrations over phi1, F2, and F3: A probability for the particular (d1,x2) configuration in question. Add these up for the possible x2s, and you get a probability for the given d1 value. Get the probs for all d1 values and you can compute the RMS value of that distribution, half of that is the sigma_y(1) Allan dev.

Assuming this is a good model, the sigma_y(1) Allan deviation value should tell us, quite accurately, what the most likely standard deviation of the 1s-window oscillator frequency F is. And similarly, looking at the y values should tell us the mean frequency (although we know this drifts over time).

I think my next step is to write a short C or Matlab or Mathematica program to actually do this calculation. Only advantage of Mathematica: Maybe we can use it to evaluate the integral analytically. If not, I can do it numerically.

Journal Quest

Main task today: Working on some write-up of my recent analysis for the journal paper.

Gordon came by and I showed him the recent stuff. We looked at how to determine the time uncertainty from the Allan deviation. That would be a good thing to include in the paper.

Met with Ray, we went over the status of everything in the COSMICi system (HW at least).

Emailed Brian Stearns (contact at DeLorme) to ask him if the 0.2 ppm frequency instability we are seeing could be attributable to the TCXO used in the kit, requested part # for datasheet lookup. Also told him about the (10/20/200/400) millisecond glitches.

I'm wondering now if it is really a good assumption for us to say that quantization errors of neighboring phase readings are independent. Suppose the real OCXO frequency is close to an integer; the the quantization errors will be almost totally correlated, and will vanish from the y values. Or, if the OCXO frequency is close to a half-integer, then the quantization errors will be offset from each other by 0.5, and will tend to add up constructively in the y value, but then will tend to cancel again in the d (deviation) determination.

I'm thinking I need to make a more detailed model of the quantization error. Some parameters:

a - quantization error for first (x1) phase measurement, in cycles, in range [-0.5, +0.5). (That is, actual phase at the time of the PPS edge was the measured value minus a.)

f1 - fractional part of actual average frequency (phase accum. btw. PPS edges) during 1st second, in range [0, 1).

b - can be derived from f1 and a

f2 - fractional part of actual average frequency (phase accum. btw. PPS edges) during 2nd second, in range [0, 1). This can be derived from f1 and a model of the "real" instability, as a normal distribution.

c - can be derived from f2 and b.

Then, y1 and y2 can be derived from this, their error (relative to f1 and f2) can be derived, and the d=y2-y1 can be derived, and the error between this and the real delta of f2-f1 can be derived.

We can sweep through values of f1 and a, and plot the distribution of deviation measurement, including the effect of quantization error, as a function of the "real" underlying deviation.

This is all well and good, but I think the argument for phase offsets being pretty well independent for larger window sizes is pretty strong, otherwise I think we'd see some variation in the adjusted Allan deviation measurements depending on the window size, which we don't see except for very small window sizes. Still, doing the above more detailed calculation would be a good sanity check on our earlier calculations.

Tying Loose Ends

The rest of last week was just spent wrapping up bits of the analysis from Monday. Getting ready to start writing. Oh, and the XMath PO finally got generated.

Mathematica FAIL

* Ray & I discussed strategy w.r.t. how to proceed with FEDM board testing. We'll wait a few more days to hopefully get documentation from Sachin/Vaibhav, since they have been ill.

* Took forms to Lee hall to get Xmath past EIT approval. The secretary said, the guy who does the approvals is supposed to be in tomorrow.

* Since the result from Mathematica was wrong, I manually constructed a symmetrized version of the distribution of frequency drift values we'd expect to see as a result of quantization error; & updated Monday's notes accordingly. Today planning to numerically calculate the distribution in C to make sure the curve is correct.

* Checked our NetBeans installation to see if it had C/C++ support installed. No; apparently I only installed it with Java support before. Uninstalling & reinstalling to add support for other languages. In the meantime, maybe I should just use Python...

* Had trouble registering NetBeans with Oracle, so Python won the day. Here is my Python code:

def dist(d):
if d > 100:
return 0.25*(2-(d/100))**2;
elif d < -100:
return 0.25*(2+(d/100))**2;
else:
return 0.25*(2-(d/100)**2);

count = {}

for d in range(-200, 200):
count[d] = 0

for x1 in range(-50,50):
for x2 in range(-50,50):
for x3 in range(-50,50):
y1 = x2 - x1
y2 = x3 - x2
d = y2 - y1
count[d] = count[d] + 1

tot = 0
totd = 0
for d in range(-200,200):
print(d, count[d], dist(d), sep='\t')
tot = tot + count[d]
totd = totd + dist(d)

print("TOT", tot, totd)

The two columns, when normalized, are identical, except for tiny discrepancies which can be chalked up to the imprecision that's inherent to the numerical integration. So, this confirms my previous belief that the solution from Mathematica was correct except for the gap.

Just call me a quant...

In the car on Friday I was thinking about analyzing the quantization error a little more precisely, so that we can more accurately subtract it out from other influences on the Allan deviation to get a clearer idea of how much of the deviation is due to actual frequency instability.

The nice thing about quantization error is that it can be well modeled by a simple uniform distribution over the interval [-50 ns, +50 ns], since 100 ns is the clock period of our OCXO, and we are doing a simple single-edge triggered input capture on the PPS signal from the GPS.

Actually, we could improve this to [-25 ns, +25 ns] by using a dual edge-triggered flip-flop: I found a VHDL gate-level structural design for this online the other day. I tested it just now and it works. Let's put that off until later, however...

For now, we can express the distribution over the value of a single y(t_1, t_2)=x(t_2)-x(t_1) between two adjacent times t_1 & t_2 (1 second apart) by integrating a delta function as follows (with x_1=x(t_1), x_2=x(t_2)):

We should be able to solve this in Mathematica. Actually, I know what the answer is already from visualizing it geometrically; it's a sawtooth distribution. Mathematica confirms it (here using the scale 1=50ns):


But really, what we are interested in is the distribution over values of the frequency drift between adjacent windows, d(t_1) = y(t_2)-y(t_1) = x(t_1) - 2x(t_2) + x(t_3). So then this is the triple integral:

p(d) = \int_{x_1=-50 ns}^{50 ns} \int_{x_2=-50 ns}^{50 ns} \int_{x_3=-50 ns}^{50 ns} \delta(d - (x_1-2x_2+x_3)) dx_3 dx_2 dx_1, or (with y=d, a=x_1, b=x_2, c=x_3, and letting 1/2=50ns):
I found this too hard to solve exactly in my head, but I got (a trial version of) Mathematica to evaluate it for me using this expression:

Simplify[Integrate[DiracDelta[d - (a - 2*b + c)], {a, -1/2, 1/2}, {b, -1/2, 1/2}, {c, -1/2, 1/2}]]

The answer should be a piecewise quadratic function in four parts, restricted to the range [-2,+2]. Unfortunately, Mathematica does not solve it properly - it has a missing piece:

Of course, the correct answer should be symmetrical around d=0. We could fill in the missing piece in [-2,-1] by mirroring the piece in [1,2]. Doing that, we get the answer:


The middle two cases can also be expressed as (d-2)^2/4, (d+2)^2/4. We can check the normalization by integrating it over [-2,2], and we get the expected result (1), so this is a real probability distribution.

Even with just three points, this is already starting to take on the character of a normal distribution, as would be expected from adding/subtracting independent random variables.

Anyway, we can model each experimental "d" sample as a sum of a sample d from the above distribution, plus a d value obtained from a normal-distribution model of the frequency deviations in each second.

However, we note that for larger window sizes, the magnitude of the quantization error stays the same (since no matter how large the windows get, the only thing that matters for the quantization error is still just the quantization error for the 3 points at the window boundaries), so it is effectively scaled down by a factor of the window size. Meanwhile, the magnitude of the inherent deviation of the frequency means (as an effect of sample size) scales down with the square root of window size. So, for larger window sizes the quantization error becomes insignificant.

Still, let's see if we can estimate its contribution. If we treat the above quadratic distribution as if it were a statistical population, we can find its variance by integrating y^2 times this distribution for y in [-2,2], and that result is just 1/2. So the standard deviation of this distribution is just 1/sqrt(2), or 0.7071, meaning (with 1=100ns) 70.71 ns.

Using the quadrature rule for the sum of the "real" frequency deviation d(real) and the quantization error d(quant), we have that the measured d is:

d^2 = d(real)^2 + d(quant)^2

And since d(quant) = 0.7071e-7, d(quant)^2 is 0.000,000,000,000,005 = 0.5e-14 sec^2 (exactly, it turns out).

So, we can write:

d(real)^2 = d^2 - 0.5e-14 sec^2
d(real) = sqrt(d^2 - 0.5e-14)

I put this calculation in my spreadsheet, and found that even in the worst case (1-second windows), the quantization error only accounts for a mere 4.2% of the deviation measured! And the proportion goes downwards from there inversely with the window size. What is happening?

Let's look at the worst case (window size 1). The "corrected" (inferred) frequency deviation sigma is about 167 ns. The standard deviation of the frequency deviation caused by quantization error is (as mentioned above) 70.7 ns, divided by sqrt(2) because of the factor of 1/2 in the Allan deviation expression, gives 50 ns. Putting these together in quadrature, the expected overall deviation is sqrt(166^2 + 50^2) = 174 ns, which is the measured deviation. In other words, because of the way errors combine in quadrature, even though the expected magnitude of the quantization error, considered by itself, is 50/174 = 29% of the measured deviation, it cannot by itself account for more than about 4% of the measured deviation.