Unfortunately I don't have access to a frequency standard which would be accurate enough for proper comparisons. What I can do, however, is look at the Allen deviation referenced to the GPS PPS.
I measured the OCXO as free running twice and then running with the PLL controller. The PLL controller was configured to estimate the phase error every 8 seconds. The control law was then set to a time constant of 48 cycles - leading to a total time constant of around 384 seconds for the PLL.
Modified Allan deviation of the OCXO as free running (2 measurements) and under PLL control
In the short time scales the Allan deviation is seen to be approximately the same between the free running oscillator measurement and the PLL controlled oscillator measurement. This is desirable, as in the short time scales it is expected that the free running OCXO has smaller phase noise than the GPS PPS. Since the stability here is measured against the GPS PPS, the true short time scale behavior cannot be observed as it is corrupted by the noise of the GPS PPS itself.
On the other hand, going to larger time scales, we know that the GPS PPS noise becomes very low. In those time scales we see that the PLL controlled oscillator Allan deviation keeps decreasing as it follows the GPS PPS ever closer, while running free, the OCXO would start exhibiting an increase in the deviation.
For the long time scale it's quite easy to get to a low deviation with just about any PLL implementation locked to the GPS PPS. It seems to be more difficult to not amplify the noise too much in the short time-scale. Even after several iterations and a lot of tuning, there is some noise gain left. I am however happy with this for now. Perhaps in the future I'll revisit this by measuring it against an atomic reference, but for my home lab needs, this is most likely more than good enough.
Phase noise characterization
At
the time of completing the build I happened to have access to a lab
with a Rohde & Schwarz FSW26 spectrum analyzer. I used it to measure
the phase noise behavior of the device.
Phase noise from 10 Hz to 10 MHz
The OCXO datasheet gives the typical phase noise values as
-80 dBc/Hz @ 1 Hz
-120 dBc/Hz @ 10 Hz
-140 dBc/Hz @ 100 Hz
-145 dBc/Hz @ 1 kHz
-150 dBc/Hz @ 10 kHz
I
did not measure down to 1 Hz or below due to the sheer amount of time
it would have taken, and I only had day time access to said lab. It would seem likely, however, that the unit is
better than -80 dBc/Hz at 1 Hz. Noise density of -120 dBc/Hz at 10 Hz
appears to be met. However, for the higher frequencies the noise
performance is clearly not reached. This discrepancy can either be
because of the OCXO itself or due to additive noise from the 74HC04s
used as clock buffers. The result is still pretty good, so I don't mind.
The
period RMS jitter can be calculated from the phase noise through
correlating the phase noise at time t with the phase noise at t +
period. This can be integrated from the phase noise spectrum. We get
50 ns: RMS jitter = 11.48 ps
100 ns: RMS jitter = 14.88 ps
To
sanity check the calculations, I used a Rohde & Schwarz RTP164
oscilloscope - also present in the lab - to directly measure the period
jitter at different time offsets
50 ns: RMS jitter = 15.53 ps
100 ns: RMS jitter = 21.05 ps
10 ms: RMS jitter = 18.36 ps
100 ms: RMS jitter = 17.88 ps
Jitter
performance is thus quite uniform over a large correlation time range
and matches the order of magnitude computed from the phase noise.
While
I wasn't targeting a low jitter clock source, I think the jitter
performance is easily good enough for anything I could consider doing at
my home lab currently. I'm extremely pleased with the results.
Just a quick post about recovering a signal, which has been corrupted by quantization error. This is something I discovered while working on my GPSDO project, but has general applicability. This all might be common knowledge and trivial, but I enjoyed coming up with my own solution. Documenting it here so I don't forget about it.
Introduction
Consider a continuous observable, which is sampled by a low resolution process. By low resolution, I mean that there is a limited number of steps that can be differentiated. In a noiseless ideal case the response would look something like shown in figure 1.
Figure 1. Ideal quantization of continuous observable.
However, consider then that the observable is corrupted by additive noise. This causes the sampled signal to also exhibit noise. The result might look as shown in figure 2.
Figure 2. Effect of additive noise on sampling.
This noise in the sampled signal can be used as a way to increase resolution, where several samples are averaged to produce a single sample. Noise is sometimes deliberately injected into measurement signals to take advantage of this. With enough samples averaged, not only is the additive noise reduced, so is also the quantization error. This is illustrated in figure 3.
Figure 3. Result of averaging noisy sampled data.
Now, the example above exhibits a large amount of additive noise. Noise standard deviation was taken as 40% of quantization resolution. However, if the additive noise is small compared to the quantization resolution, this breaks down. See figure 4, where the noise standard deviation is just 15% of quantization resolution.
Figure 4. Effect of small additive noise.
Suddenly, the quantization error is no longer removed through averaging and a step-like pattern emerges.
An important thing of note is that there are again some continuous values for which there is only a single sampled value. It seems rather clear that in general no further resolution can be extracted for those values. However, for the following discussion it's more interesting to observe the values that do still exhibit two distinct sampled values. This is also the situation I have in my GPSDO project that motivated all this.
In the quantization transition regions averaging clearly gives more resolution, but the averaged values do not match the true values. The obvious question is whether this can be improved.
Bayes saves the day
Initially I took a Bayesian approach to the problem. When presented with the quantized value \(q\), what can be inferred of the continuous value \(c\)? This is answered by Bayes' theorem. \[ p(c | q) = \frac{p(q | c) p(c)}{p(q)}. \]
That is, the probability of the continuous variable having value \(c\) after we've observed quantized value \(q\) is proportional to the prior probability of observing \(c\) times the probability of observing quantized value \(q\) given that the continuous variable has value \(c\), which is also called the likelihood. The probability of observing the quantized value \(q\) in the denominator can be thought of just a normalizing factor.
The formula expresses how we update our knowledge of the continuous variable, when presented with a measurement. The posterior distribution is proportional to the prior times the likelihood.
Before we can formulate the likelihood, we must properly define the additive noise and the quantization. To keep things simple, let's take almost the simplest quantization possible: the round function. Also let's take the noise as normally distributed and zero mean \(n \sim N(0,\sigma) \).
The likelihood is the probability of sampling a value, given that the continuous variable has some value. With our choice of quantization, we observe quantized value \(q\) if the continuous value plus noise \(c + n\) is in \([q-0.5, q+0.5)\). The probability can be evaluated as \[ p(q|c) = p(q-0.5 \leq c+n \leq q+0.5) = p(q-0.5-c \leq n \leq q+0.5-c)\]. Which, since the noise is normally distributed, can be further evaluated as \[ p(q|c) = F_n(q+0.5-c) - F_n(q-0.5-c) = \frac{1}{2} \mathrm{erf}(\frac{q+0.5-c}{\sqrt{2}\sigma}) - \frac{1}{2} \mathrm{erf}(\frac{q-0.5-c}{\sqrt{2}\sigma}). \]
To gain more intuition, let's look at what the likelihoods look like. Let's take the noise standard deviation as \(\sigma = 0.1\). This is shown in figure 5.
Figure 5. Likelihood functions for observing quantized value 0 and 1.
Near the continuous values 0 and 1 there is an almost 100% probability for observing quantized values 0 and 1, which causes the stairs we saw in figure 4. The significant overlap area of the probabilities is only a small fraction of the range, which leads to the average not matching the true value also seen in figure 4.
Then, say that we have observed a pattern 1, 0, 0, 0, 1, 0, 0, 0. What can we say about the continuous variable? Taking the prior distribution as uniform (i.e. no knowledge of the continuous variable), the Bayesian inference progresses as shown in figure 6.
This approach is very powerful, but also very heavy. And as an additional complexity, if the system is dynamic, we need to add uncertainty to the distribution between inference steps, to allow the inferred continuous variable to evolve. For the usual applications, this likely adds too much complexity to be practical.
This exercise, however, gave me insight in why the simple averages of the sampled values don't give the correct result, as well as how the likelihoods look like.
A practical solution
In practice, the sampled values alternate at most between 2 distinct values. While it's not strictly impossible to observe 3 or more values, for the situations of interest here, the probability is just extremely low. This allows us to simplify things by allowing only two values. The likelihood of observing a zero is then simplified as \[ p(q=0|c) = F_n(0.5-c) = \frac{1}{2} + \frac{1}{2}
\mathrm{erf}(\frac{0.5-c}{\sqrt{2}\sigma}), \]while the likelihood of observing a 1 is the complement of this. These likelihood functions are visualized in figure 7.
Figure 7. Likelihoods of the two observable values.
Since the two cases are complements of each other, inferring the continuous variable becomes much easier: probability of observing, say, the value 1, can be directly mapped to the value of the continuous variable. The probability, on the other hand (for the values 0 and 1), can be simply estimated as the average of observed samples. Figure 8 shows the mapping from probability to value.
Figure 8. Continuous variable inference from probability of observed value.
While it appears in figure 8 that a probability of 0 gets mapped the continuous variable value 0 and that the probability of 1 gets mapped to the continuous variable value 1, this is not quite the case. This is because there are still very tiny probabilities unaccounted for in figure 8. However, since these are practically impossible, we can round the remaining probability to zero.
Since estimating the probability of observing a 1 is the same thing as computing the average of ones and zeros, we can extend the concept of averaging the values to the other values as well. Turns out this way we can define a single curve that takes the average of the observations and outputs an estimate of the continuous variable. See figure 9.
Figure 9. Continuous variable estimation from average of observations.
This approach is very lightweight to compute, and allows handling dynamic situations by choosing the length of time that the average is computed over.
Let's revisit the situation of figure 4 using this approach. The results are presented in figure 10.
After the good results I got from my prototype GPSDO, I wanted to go a bit further than the mess of wires that was the prototype. Besides achieving a better build quality, I needed more than just one 10 MHz output to support a few key instruments in my lab - and a PPS output would also be nice to have. I also wanted to improve the temperature dependency of the device by insulating it better, and by keeping the control voltage reference at constant temperature.
The device was built in April - June of 2024, but I only now got to documenting it.
Overview
The device is composed of two parts: the main unit and the GNSS receiver unit. They are connected to each other through standard twisted pair Ethernet cabling, which allows placing the receiver unit to a convenient location, which may be far away from the main unit.
GNSS receiver unit on the left, main unit on the right
Reverse sides
GNSS receiver unit
The GNSS receiver unit gets power through one of the four pairs of the cable. A suitable voltage for the receiver is then regulated locally from the provided voltage. The PPS output of the receiver is fed into a RS422 line driver and transmitted back to the main unit over another pair of the cable.
The receiver serial TX and RX are also connected to a RS422 transceiver to allow data communication with the main unit. Thus, all 4 pairs of the cable are used.
Cover dome, board carrier and pole mount are held together with 4 M3 screws
Close up on board carrier. RJ45 is on a separate board and connects via ribbon cable
RJ45 board is held in place with snap tabs. Receiver PCB is attached with 2 M2.5 screws and hold down tabs.
Reverse sides
Main unit
The main unit is comprised of several sub-units. These are the main board sub-unit, the OCXO board sub-unit and the LCD board sub-unit. Each sub-unit handles a separate part of the system and communicate with each other through I2C.
Main board sub-unit and LCD board sub-unit are attached to the cover plate. OCXO board sub-unit is inside block of foam insulation.
Mess of wires for the connectors. All wires connect to the main board. RJ45 connector is on a separate board, which is held in place by snap tabs.
LCD board sub-unit (top) and main board sub-unit (bottom)
Block of foam insulation with OCXO board sub unit inside removed from the main unit enclosure
LCD board sub-unit
The LCD board sub-unit handles displaying information to the user through a HD44780 type liquid crystal display. It acts as an I2C master and pulls information to display from the main board sub-unit as well as from the OCXO board sub-unit. It is built around a CH32V003 MCU.
Close up on the LCD board sub-unit with main board sub-unit
OXCO board sub-unit
The OCXO board sub-unit handles the 10 MHz frequency generation and its control voltage generation, and thus carries the OSC5A2B02 OCXO itself. It is housed inside a thermally insulating foam block to reduce the temperature dependency of the 10 MHz output.
The OCXO adjustment is about 1 ppb per millivolt, so the control voltage needs to be quite stable. This stability is provided by a TL431C voltage reference, which is at close thermal contact with the OCXO and thus also at a constant temperature. A 16 bit DAC is then implemented using a PWM signal from a CH32V003, followed by heavy filtering to try to keep the phase noise low.
The OCXO board sub-unit also features a DS18B20 temperature sensor, which is at close thermal contact with the OCXO. This is for now only used for monitoring, but could perhaps be used for some additional temperature compensation in the future.
The 10 MHz output from the OCXO is not used by any of the circuitry on the sub-unit itself, but is instead passed to the main board sub-unit.
The OCXO board sub-unit acts as an I2C slave and allows setting the OCXO control voltage based on a control word, which it receives from the main board sub-unit. Additionally, the LCD board sub-unit queries it for temperature information.
Block of foam insulation with OCXO board sub-unit inside
Foam cover removed
OCXO board sub-unit removed from block of insulation
Main board sub-unit
The main board sub-unit performs three functions: GNSS data interfacing, OCXO PLL control and clock buffering.
Main board
Reverse side
GNSS data interface
The GNSS data interface has a RS422 transceiver to allow serial data communication with the GNSS receiver. This is used to configure the receiver, as well as for receiving a time stamp for each PPS pulse for absolute phase control. The GNSS data interface is built around a CH32V003 and acts as both an I2C master and an I2C slave. As a master, it actively pushes validity and time stamp information to the OCXO phase lock loop system. As a slave it allows the LCD board sub-unit to fetch the validity, time stamp, positioning and other GNSS information.
PFD controller
The phase-frequency detector controller is the heart of the system, and is built around a CH32V003. The MCU is clocked from the OCXO board sub-unit 10 MHz output. This 10 MHz is internally clock doubled for a 20 MHz cycle rate at the MCU core.
Via an RS422 receiver, the MCU gets the PPS signal from the GNSS unit. This allows the MCU to count the number of 10 MHz cycles between PPS pulses at half-cycle (i.e. 50 nanosecond) resolution. This, averaged over a long period, represents the true frequency of the OCXO output.
At each PPS, the accumulated phase count of the OCXO is compared to the ideal phase (i.e. the phase derived from the time stamp). The difference of these represents the phase error. The phase error is then propagated through the PLL software controller to compute an OCXO control word.
The OCXO PFD controller MCU acts both as an I2C master and as an I2C slave. As a master, it transmits the OCXO control word to the OCXO board sub-unit. As a slave it allows the LCD board sub-unit to fetch status information regarding the PLL lock.
Clock buffers
The clock buffers allow isolating the external users of the 10 MHz clock signal from the clock generation. Thus the PLL keeps running even if some clock outputs were e.g. accidentally shorted or otherwise abused. There are a total of four 10 MHz outputs and one PPS output.
Each clock buffer is made from a 74HC04 hex inverter. Five of the inverters are paralleled for the output stage. This is done to provide a strong drive for a 50 ohm output. The sixth inverter is used as a buffer between the OCXO signal and the five paralleled inverters. This reduces the capacitance which the OCXO needs to drive as it now only sees the capacitance of one CMOS input per clock buffer instead of six.
Characterization
I think I've now spent more time measuring and tuning the device than actually building and designing it. In fact, I'm not yet even using it as a time base for any of my lab instruments.
OCXO gain
The OCXO is modeled as having a linear response to its control voltage (see previous post), while the control voltage is assumed to be linear to the control word. The OCXO frequency change per control word change is here called the OCXO gain. In order to tune the controller for the best response, the OCXO gain needs to be measured at the operating point.
OCXO gain measurement
The control word given to the OCXO board sub-assembly was varied between two values, one of which was 1000 counts above the operating point and one which was 1000 counts below the operating point. The frequency of the OCXO was then measured against the GNSS PPS signal. Note that since the PFD operates with a clock-doubled OCXO input, the gain is defined with respect to the clock-doubled 20 MHz frequency and not the 10 MHz output.
The clock-doubled OCXO output changed frequency by 0.458 Hz with a control word change of 2000, thus giving the OCXO gain an approximate value of 229 micro-Hz per count.
OCXO PLL time-scale
A critical parameter to choose is the bandwidth of the PLL. This is because the GNSS PPS is noisy at short time scales, while the OCXO is unstable at long time scales. To gain understanding on the matter, I collected phase error data from the OCXO while set as free running.
As comparison points, I found some stability data published for the OSC5A2B02 as well as for bare GPS PPS signals. These came from the blogs PA1EJO and www.febo.com. These sources used rubidium standards as their reference clocks, while I could only compare against the GNSS PPS.
Phase error of the OCXO against GNSS PPS
Low frequency phase noise on the 10 MHz output
Modified Allan deviations. My measurement and other published data.
Looking especially at the modified Allan deviation graphs, my data follows the [PA1EJO] GPS PPS data quite well at the low frequencies. This indicates that the deviation at those frequencies comes from the PPS and not from the OCXO. On the other hand, at slightly higher frequencies my data follows the [www.febo.com] OSC5A2B02 data, indicating that here the OCXO is causing the deviation. Though my OCXO in this measurement appears to be more stable than the unit in the external data. This could be due to additional thermal insulation my unit has.
The OCXO has apparent linear frequency shift, which at such long time scale could be easily corrected by the PLL. To still get more understanding of the stability, I computed the Hadamard deviation as well. This is insensitive to linear frequency drift, and the [PA1EJO] GPS PPS data has this metric also published.
Hadamard deviation comparison
In the Hadamard sense, the OCXO is remarkably stable. There is a slight change in the slope at around 10 seconds time scale, but it isn't anything to really worry about.
It looks that it's only important to select the time scale short enough, that the linear frequency drift is compensated, but still long enough that phase error measurements (through low-pass filtering) become accurate. This is probably a pretty wide time scale range. Mostly it's trying to get the scale as short as possible, while keeping the control stable.
To explain the last part a bit. The hardware can directly observe the phase error at only 50 ns resolution. This by itself is not nearly good enough for proper control. However, due to jitter in the PPS, the phase error is actually observed alternating between two values. Low-pass filtering this alternating raw error gives much improved resolution. As an additional trick, the controller can provoke jitter in the measurement by deliberately controlling the OCXO phase to lie at a midpoint between two values.
Conclusion
I'm still in the progress of measuring the response of parameter choices, as well as making tweaks to the controller algorithm itself. There will be a part 3, in which I'll try to get some performance measurements. I'll also try to get the electronics and code published as well.
In a previous post I described some of my experiments with building a GPS disciplined oscillator. In my experiments, I lock a voltage controlled oscillator to the PPS signal of a GPS receiver. The locking is achieved by a control loop, which observes the phase error between GPS and the VCO and produces a control signal to adjust the frequency of the VCO.
The GPS PPS signal is very noisy at higher frequencies (above say 1 milli-Hz), so the controller must have a low bandwidth to reject the higher frequency noise. Also, since I don't want to spend a lot of time tuning the controller, I want the bandwidth to be the only parameter.
VCO phase error model
The VCO has a control voltage input \(u\), which controls its frequency \(f_\text{vco}\). This relationship is well approximated as linear around an operating point. Thus near our target frequency \(f_\text{target}\) we can model the VCO frequency as
\[f_\text{vco}(t) = f_\text{target} + g (u(t) - u_o),\]
in which \(g\) is the VCO gain and \(u_o\) is the needed control value to attain the target frequency. The frequency error is thus given by
\[ f(t) = f_\text{target} - f_\text{ocxo}(t) = -g (u(t) - u_o) \]
As the goal is to achieve phase lock, the quantity of interest is actually the phase error. We take the unit of phase to be a full cycle. This makes the phase error \(e\) simply the integral of the frequency error over time.
\[e'(t) = - g (u(t) - u_o) \]
As the controller will be implemented in the digital domain, \(u\) will be piecewise constant. This allows us to discretize the model as the recurrence
\[ e_{n+1} = e_n - \Delta t g (u_n - u_o),\]
in which \(\Delta t\) is the time between the discrete changes of \(u\), \(u_n\) is the value of \(u\) at \(t \in [n \Delta t, (n+1) \Delta t) \) and \(e_n = e(n \Delta t)\).
Controller
The simplest form of controller which can control the phase error to zero is a PI controller. Here the control is determined as a linear combination of the phase error \(e\) and its integral \(i\) as
\[ u_n = P e_n + I i_n, \]
in which \(P\) and \(I\) are tuning parameters of the controller.
Unfortunately for us, our phase error measurement is very noisy at short time scales. With the simple PI controller the noise would couple straight to the control output through the P term. We thus need to consider something else. One could obviously use a longer interval between the measurements, but that introduces its own set of problems. Our solution is to instead apply a simple first order low-pass filter to the phase error measurement and then combine that with a PI controller.
Since the filter and the controller are implemented in the digital domain, let's transition now to an entirely discretized domain and ignore the continuous domain altogether. In this context, we take our low-pass filter as
\[ \hat{e}_{n+1} = (1-\alpha) \hat{e}_n + \alpha e_n, \]
in which \(\hat{e}\) is the filtered phase error and \(0 < \alpha < 1\) is a parameter defining the filter bandwidth.
The integral of the filtered phase error, similarly, is considered only in the discretized sense, and is given as
\[ \hat{i}_{n+1} = \hat{i}_n + \hat{e}_n \]
The control law of the controller is then
\[ u_n = P\hat{e}_n + I\hat{i}_n, \]
in which \(P\) and \(I\) are tuning parameters to set the behavior of the controller.
Combining the VCO phase error model with the controller gives the model of the entire system as
\[ \begin{equation}
\begin{pmatrix} e_{n+1} \\ \hat{e}_{n+1} \\ \hat{i}_{n+1} \end{pmatrix} = \begin{pmatrix} 1 & -\Delta t g P & -\Delta t g I \\ \alpha & 1-\alpha & 0 \\ 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} e_n \\ \hat{e}_n \\ \hat{i}_n \end{pmatrix} + \begin{pmatrix} \Delta t g u_o \\ 0 \\ 0 \end{pmatrix}
\end{equation} \]
For the controller to work, we want the errors to converge. The dynamics of the convergence are determined by the eigenvalues of the matrix. For simplicity, we'll choose all three eigenvalues as \(r\), where \(0 < r < 1\).
The characteristic polynomial of the matrix is
\[ \lambda^3 + (\alpha - 3) \lambda^2 + (\Delta t g P \alpha - 2\alpha + 3) \lambda + \Delta t g I \alpha - \Delta t g P \alpha + \alpha - 1. \]
On the other hand, to have the desired eigenvalues, we want the characteristic polynomial to be
\[ \lambda^3 - 3r\lambda^2 + 3r^2 \lambda - r^3. \]
Matching the coefficients, we get the equations
\[ \begin{align} \alpha - 3 &= -3r \\ \Delta t g P \alpha - 2\alpha + 3 &= 3r^2 \\ \Delta t g I \alpha - \Delta t g P \alpha + \alpha - 1 &= -r^3 \end{align} \]
Solving those equations yield
\[ \begin{align} \alpha &= 3(1-r) \\ P &= \frac{1-r}{\Delta t g} \\ I &= \frac{(1-r)^2}{3 \Delta t g} \end{align} \]
The resulting matrix is unfortunately non-diagonalizable, making analyzing the resulting dynamics a bit tedious. The controller however appears to be well-behaved and produces convergence without too much overshoot.
The following figures show simulated responses from a controller with \( r = 0.001 \).
Impulse response of the controller against phase error change
Impulse response of the controller against control offset change
The CH32V003 reference manual does not explain some key information for using the I2C peripheral. The issue is mainly with the I2C1_CKCFGR register, but there is also confusion around the FREQ field of the I2C1_CTLR2 register.
Digging a bit, it turns out that WCH appears to be using the same I2C peripheral IP as Gigadevices and Puya use in many of their microcontrollers. Both of the other vendors have a bit more documentation in their manuals.
The explanation of the I2C1_CKCFGR register turns out to be fairly simple, and the values used by the I2C examples going around the net appear to be correct. Its documentation in the reference manual says the following
However, what is left unsaid is the following
Mode
Duty
T_high clock cycles
T_low clock cycles
F/S=0
Either
CCR
CCR
F/S=1
DUTY=0
CCR
2*CCR
F/S=1
DUTY=1
9*CCR
16*CCR
The bitrate is just BR = F_APB1 / (T_high + T_low), and thus the magical factors of 2, 3 and 25 used in the examples are explained. This also explains where the quoted duty cycles of 33% and 36% come from.
This still leaves the FREQ field in I2C1_CTLR2 to be explained. It is described as
The Gigadevices and Puya documentation agree with WCH's documentation on this: you're expected to program the field with a value, which is the peripheral clock frequency given in integer megahertz.
All examples for the CH32V003 I've seen on the net, however, program it as F_APB1 / 2000000, which is half of the correct value. Also, many of the examples document the value 2000000 as the frequency of the peripheral, which is not correct at all. WCH's documentation is worded in less obvious terms than the other vendors, which may be the source of the confusion.
But what does this register do? If it's possible to write only half of the correct value to the field and still have the peripheral seemingly work, it can't be very critical. This is only covered in the Puya documentation, and still with not enough detail. The English translation states:
This register must be configured with the value of the APB clock frequency to generate data setup and hold times that are compatible with the I2C protocol.
Not sure why the timings given in the I2C1_CKCFGR are not enough to guarantee proper setup and hold, but apparently there is additional signal conditioning which this field affects.
The best (overall) frequency reference I have is in my Racal-Dana 1992 frequency counter, which has the ovenized oscillator option 04A. It's great, but calibrating it with the equipment I have is an annoying task. I calibrate it against GPS. However, as GPS time is stable enough only in the long term, it requires measuring a lot of GPS PPS samples. This takes a long time, and even longer if I want to adjust the time base. This got me thinking if I could automate the process, and after some thought I realized the best approach would be a GPS disciplined oscillator.
While there are commercial GPSDO devices for sale, I can't justify the cost for buying one. I could however design and build one. Turns out I can justify much higher cost for such a project, as there's a a lot of benefit in the learning experience.
My goal would be to produce 10MHz +- 1ppb. I also want to keep the phase noise low, but that's somewhat secondary and without quantitative requirements.
The GPS modules I've used are really old Fastrax parts. They are specified to 50ns RMS jitter on the PPS output. Not great, but workable. This pushes the integration time to around 1 hour to get down to 1 ns jitter. I would not trust the PPS pulse below a 1 hour window anyway (more on that later). Anyway, this puts a requirement on the stability of the oscillator: it would need to be stable (better than +-1ppb) for much longer than the GPS integration period - say 10 hours. This puts me well into the OCXO domain.
DIY OCXO - fail
I happen to have some old VCTCXOs (voltage controlled temperature compensated crystal oscillator), which I found thrown out years ago. They seemed high quality and expensive (from Rakon), so I always wanted to find a use for them. I first thought to build temperature control around one of those and convert it to an OCXO. Using a DS18B20 and PWM controlling a heater I managed to get the temperature to be controlled well within 0.1 degrees Celsius. However, I could not get the crystal oscillator to behave though. I kept the control voltage grounded to make sure any noise on the control voltage wasn't causing the issue. In the end it may have been due to the supply voltage sensitivity of the part, which is not too good at +-300 ppb for +-5% supply voltage. This would require about 1 mV stability for the supply voltage for my application.
Rakon VCTCXO frequency and temperature, moving average over 5 minutes Measured with Racal-Dana 1992 and DS18B20
Looking around the internets, I found that all the cool cats are playing with OSC5A2B02 OCXOs, which are available on Aliexpress for very cheap. So I ordered some. At about 3€ a piece it didn't seem like a too big investment.
OSC5A2B02 testing
The OSC5A2B02 is much less sensitive to supply variations (+-2ppb for +-5% supply voltage), so just about any regulator suffices. What is important to note however is the control voltage sensitivity, which is about 1000 ppb per volt, or equivalently 1ppb per millivolt. For the initial tests, I simply grounded the control voltage to eliminate it's contribution.
OSC5A2B02 frequency (insulated in PE foam), moving average over 5 minutes Measured with Racal-Dana 1992
The variation was much larger than the datasheet of the oscillator promised. This is when I realized, that I actually didn't have a clear idea of what the stability of my frequency counter was, and that I could just be seeing it's variation.
To get a better idea of the frequency counter variation, and thus the true stability of the OCXO, I connected a GPS PPS on the B channel of the frequency counter.
OSC5A2B02 frequency (top) and GPS PPS (bottom), moving average over 20 minutes Measured with Racal-Dana 1992
Bingo! Turns out most of the variation observed in the OSC5A2B02 frequency is actually due to variation of my frequency counter! Who would have guessed that a 3€ Aliexpress OCXO today is so much better than an instrument worth thousands in the 1980s.
Normalizing the OCXO frequency with the PPS frequency (after heavy filtering), the OCXO appears to be within +-1ppb over a day as long as there is no control voltage variation.
Phase noise and jitter
The OCXO datasheet gives some spot values for the phase noise:
-80 dBc @ 1 Hz
-120 dBc @ 10 Hz
-140 dBc @ 100 Hz
-145 dBc @ 1 kHz
-150 dBc @ 10 kHz
Integrating over this range gives ~0.3 mrad RMS of jitter, or about 4.8 ps.
Control voltage generation
As said earlier, the control voltage is about 1 ppb per millivolt over a range of 4 volts total. Using 16 bit discretization for the range leads to about 0.06 ppb resolution. To get a ballpark figure of how least significant bit transitions affect the phase noise, let's consider the case of modulating the 10 MHz carrier by +-0.03 ppb ( = +-0.3 mHz) at 1 Hz frequency (assume single tone sinusoidal modulation for simplicity). This is given by
\[ y(t) = \cos(2 \pi f_c t + A \sin(2 \pi f_m t)) \]
In which \( A = \frac{\Delta f}{f_m} = \frac{0.3 \text{mHz}}{1 \text{Hz}} = 3 \cdot 10^{-4} \).
Using the sum of angles identity for cosine, we get
Since \( A << 1 \) we can approximate \( \cos(A \sin(2 \pi f_m t)) \approx 1 \) and \( \sin(A \sin(2 \pi f_m t)) \approx A \sin(2 \pi f_m t) \). This gives
Thus considering the single sideband phase noise, we see that the modulating frequency is attenuated by \( \frac{A}{2} \) with respect to the carrier, which is -76.5 dBc with our numbers - or an additive 3.3 ps. So not great, but not terrible - especially if we minimize the occurrence of transitions by adding some hysteresis. Anyway, we shouldn't use any less than 16 bits of resolution.
The next issue is how to implement a 16 bit DAC cheaply. Here a pretty obvious candidate is to use PWM. Proper DAC chips with 16 bits cost close to 10 euros, while PWM and a lot of filtering can be achieved with less than a euro. Much less than 1 Hz of bandwidth is perfectly fine, so the question is just how much filtering is needed. Assume the PWM repetition period is 400 Hz and amplitude is 4V. Take then the attenuation of the filter as \(G\). Recalling that the sensitivity of the OCXO is 1 ppm per volt, the power of such a modulation relative to the carrier is given by \( 0.05 G \) (Later edit. Reading back this text, I was as confused as you are about where I got this the number. The factor comes from the above calculation, where we see that modulating the phase with amplitude A is seen as a spur of relative amplitude A/2 at an offset frequency. A 4V control modulation causes a 4 ppm change in the frequency, i.e. 40 Hz. The phase modulation amplitude is thus A = 40 Hz / 400 Hz = 0.1. Half of this is seen in the output, thus 0.05 or about -26 dBc). To push the additive jitter down to the same scale as the intrinsic jitter of the oscillator, the attenuation needs to be at least 44 dB. This should be easily achievable with simple RC filters.
The main issue with the control voltage is thus the stability requirement. The control voltage reference should remain within 1 mV over several hours. This could be combated with a high stability voltage reference. Problem is that those are expensive, and I would like to keep everything as cheap as possible. As we already have a temperature control loop (in the OCXO itself that is), we might as well use that to keep a voltage reference at constant temperature also. Turns out a TL431C has a typical stability of 4 mV over the entire temperature range, and hopefully better than 1 mV when kept at constant voltage. Also turns out TL431s are really cheap.
Microcontroller
The plan is to implement the PLL using a microcontroller. The uC would be clocked from the OCXO and it times the interval between PPS pulses. This allows determining the OCXO frequency.
The OCXO works on 5V. It would be useful if the microcontroller would also operate on a 5V supply. It also needs to have hardware input capture features, allowing precise measurement of the PPS period against the OCXO frequency. And like everything else, it must be cheap. As luck would have it, I recently ordered 50 units of CH32V003 controllers for 0.20€ per piece (including shipping). They check all the boxes for this project. The even have an internal PLL to allow doubling the 10 MHz OCXO clock for improved timing resolution.
Prototype
I implemented a very quick hack of the system on a CH32V003. It is clocked from the OSC5A2B02 with an internal PLL configured to double the frequency. This clock is used to drive a timer peripheral, which is configured for input capture from the PPS pulse. Timer resolution is increased in software from the HW provided 16 bits to 32 bits - otherwise a full PPS period could not be counted.
After each captured PPS pulse, a simple validation is performed to try to ignore erroneous pulses. A simple PI controller then controls to minimize the frequency error (frequency lock loop). The controller is just parametrized by the control bandwidth and was designed to have critically damped dynamics. Controlling just the frequency still leaves an uncontrolled phase error. I'll be looking into that too, but so far it is of no concern.
The DAC is implemented with a 16 bit PWM running at 305 Hz and filtered with a ~0.2 Hz (-3 dB) first order lowpass filter. The filtering leaves a lot to desire, and this will have to be improved for the final product. The voltage reference is provided by a jellybean TL431C, which is not temperature controlled.
Crude schematic of the prototype voltage control
Operating the control voltage as shown in the figure above allows reducing the full swing range, which increases resolution.
Experiments with the prototype have shown better than expected performance. Though I'm still lacking a data channel for the GPS data, which means that I don't know when the PPS is valid, nor can I keep absolute phase. This leads to some bad samples passing validation and causing trouble.
Closed loop controlled frequency, locked to GPS PPS GPS PPS moving average over 5 minutes
The figure above shows closed loop control of the OCXO locking on to the GPS PPS frequency. The initial condition was deliberately set about 100ppb off to see the dynamics of the control. The bandwidth of the controller was chosen as ~4 mHz for this experiment to perform the experiment at a reasonable speed. This is too wide a bandwidth for proper operation though. From the plots, we see that the dynamics are nearly critically damped, with only very minor overshoot. This is good enough for me!
75 hour long term stability experiment GPS PPS moving average over 1 hour
The next experiment was a long term stability test. Here I set the controller bandwidth to ~0.1 mHz to better reject noise in the PPS signal. During the experiment, there were a few moments at which the PPS became invalid. This causes the sudden spikes seen in the error graph. Regardless, the error remained well within +-1 ppb, which was the design goal. Also, the prototype only implements very simple control voltage filtering and the voltage reference is still just at room temperature. The design goal thus seems very much achievable, and possibly going down to +-0.25 ppb is possible with more care taken. The control voltage graph is computational, based on an assumed ideal voltage. I don't know for sure if the drift seen in it is due to the OCXO really needing to be adjusted or the voltage reference drifting and needing compensation, but my guess is on the voltage reference.
Nils Berglund, who is the author of many amazing simulation videos on Youtube, recently posted videos showing visualizations of waves on a sphere. In particular, this video caught my eye, depicting tsunami waves of the 2011 Japan earthquake travelling across the Pacific ocean.
In fact, back in 2011 I had made a very similar visualization inspired by the events, and later uploaded it to Youtube - see video below. The video was made as promotional material to get new students to major in mathematics at Aalto University. At that time I was a graduate student at the maths department working toward my PhD.
Anyway, it was as an undegrad in 2007 when I first got interested in solving the wave equation on the sphere. I had been studying the finite element method with great interest, and had programmed my own simple implementation of nodal P1 elements on triangles. I had succesfully used that to solve waves in planar domains - see video below. I thought that in order to solve the wave equation on a sphere, I would just use a mesh approximating the sphere instead. Turns out it's almost as simple as that, but there are a couple nontrivial details.
My undergrad studies included a mandatory "special assignment". These assignments were often to write a research paper on a subject - though they were published only internally at the department. I suggested to my professor that as my assignment I could study solving the wave equation on spheres using FEM.
Now, finally, I thought I'd revisit the topic, and have something generic written down on the subject a bit more publicly.
Wave equation and the finite element method
To get up to speed, let's look at the wave equation and the finite element method in a quick and non-rigorous fashion. The wave equation (in strong form) is
in which \( \ddot{u} \) is the second derivative with respect to time, \( c \) is the wave speed, which we'll just take as 1 from here on, and \( \Omega \) is the domain in which we're solving the equation at. In addition, there are initial conditions and boundary conditions that need to be described, but let's skip those for now.
FEM works with the weak formulation of the problem. To get the weak form from the strong form, both sides of the equation are multiplied by a test function \( v\) and integrated over the domain \( \Omega \). Integration by parts is applied to reduce the regularity requirement of \( u \) from 2 derivatives down to 1. The weak form is then: find \( u \in V \) such that
\[ \int_\Omega \ddot{u} v\ \mathrm{d}x + \int_\Omega \nabla u \cdot \nabla v\ \mathrm{d}x = 0 \quad \forall v \in V \]
With FEM, as with other Galerkin methods, we replace the space \(V\) with a finite dimensional approximation \(V_h\). We come up with a basis \( \{v_1,\ldots,v_N\} \) of \(V_h\), and write the finite dimensional approximate solution \( u_h \) as
\[ u_h = \sum_{i=1}^N U_i v_i \]
The problem then becomes: find \( U_i \) such that
The idea in the finite element method specifically is to define the basis functions \( v_i \) using a mesh, which fills the domain \(\Omega\). Assuming now a triangular mesh, a simple approach is to take each basis function as a piecewise first order polynomial, each piece being one triangle. Then, assign \(v_i = 1\) when evaluated at vertex \(p_i\), and \(v_i = 0\) when evaluated at any other vertex.
The integrals can be written as the sum of the separate contributions from each element
For any given \(i\) and \(j\), on the vast majority of the elements, \(v_i\) or \(v_j\) (or both) is identially zero on the element, and thus they don't contribute to the integral. Contributions come from only the triangles, that contain both vertex \(p_i\) and vertex \(p_j\).
Now, to actually compute the integrals over a triangle, a standard trick is to do a change of variables and perform the integration over a reference triangle instead. A coordinate on the reference triangle \( \hat{x} \) is easily mapped to \( x \) on the target triangle \(T\) via an affine transform \( x = A_T \hat{x} + b_T \). Let's take the triangle \(\hat{T}\) with vertices \( \hat{p}_1 = (0,0) \), \( \hat{p}_2 = (1,0) \), \( \hat{p}_3 = (0,1) \) as our reference triangle. With this choice, and the fact that affine transforms map first order polynomials back to first order polynomials, the basis functions can be expressed on the reference element in a simple way.
The 3 basis functions on the reference element are such that they get the value 1 in one vertex and 0 in the others. Written explicitly, they are
Here's where things get at bit non-standard from the usual FEM implementation. Our reference element lives in \(\mathbb{R}^2\) while the global element lives in \(\mathbb{R}^3\). This means that the matrix in our affine mapping is non-square, and we need to be a bit more careful with considering how the area scales and how the functions change under the change of variables.
Usually, with the change of variables \(x = f(\hat{x})\), we see the differential area element change as \( \mathrm{d}x = |\det J_f(\hat{x})| \), where \( J_f(\hat{x}) \) is the Jacobian of \( f(\hat{x}) \). This does not however make sense for a non-square \(J_f\) and thus something more general is needed. Turns out, that a slightly more general expression for the differential area element is \( \mathrm{d}x = \sqrt{\det J_f(\hat{x})^t J_f(\hat{x}) }\ \mathrm{d}\hat{x} \). Since \(J_f^t J_f\) is always square, this is well-defined in our case. Note, that in the case where \(J_f\) is square, this reduces to exactly the usual form. In our case, the Jacobian of the mapping is simply the matrix \(A_T\) of the affine transform, and we have \(\mathrm{d}x = \sqrt{\det A_T^t A_T}\ \mathrm{d}\hat{x} \).
The contribution of element \(T\) to \(M_{ij}\) can be computed as
While that's all for the mass matrix, for the stiffness matrix we have the additional headache of how the gradient will transform under the mapping. For basis function \(v_i\) we have \(v_i(x) = v_i(A_T \hat{x} + b_T) = \hat{v}_\hat{i}(\hat{x})\). Skipping a bit of detail for brevity, it follows from the chain rule that \( (\nabla \hat{v}_\hat{i})(\hat{x}) = A_T^t (\nabla v_i)(x) \), which by itself is enough for the typical case as \(A_T^t\) is invertible. For our case, however, we must consider additionally that the gradient \((\nabla v_i)(x)\) is actually a vector in the tangent space of \(\Omega\) at \(x\). Since the tangent space of \(\Omega\) within element \(T\) is spanned by the range of \(A_T\), there must exist some \(g\) such that \((\nabla v_i)(x) = A_T g\). Thus
\[ A_T^t A_T g = (\nabla \hat{v}_\hat{i})(\hat{x}) \]
from which it follows that
\[ g = (A_T^t A_T)^{-1} (\nabla \hat{v}_\hat{i})(\hat{x}) \]
After spatial discretization with FEM we're left with the ordinary differential equation
\[ M \ddot{U} + K U = 0 \]
Without going into detail, the problem is particularly well suited for e.g. the Crank-Nicolson method or the leapfrog integration method. These methods are stable for undamped oscillating problems and they do not produce decaying solutions. The key word for such methods is symplectic.
Crank-Nicolson
For Crank-Nicolson, we'll first convert the equation into first order form as
\[ \begin{pmatrix} I & 0 \\ 0 & M \end{pmatrix} \begin{pmatrix} \dot{U} \\ \ddot{U} \end{pmatrix} = \begin{pmatrix} 0 & I \\ -K & 0 \end{pmatrix} \begin{pmatrix} U \\ \dot{U} \end{pmatrix} \]
Denote then
\[ B = \begin{pmatrix} I & 0 \\ 0 & M \end{pmatrix} \]
\[ A = \begin{pmatrix} 0 & I \\ -K & 0 \end{pmatrix} \]
and
\[ z = \begin{pmatrix} U \\ \dot{U} \end{pmatrix} \]
to get
\[ B \dot{z} = A z \]
For the discretization, consider the derivative on the left size replaced by the approximation \( \frac{z_{n+1} - z_n}{\Delta t} \), and consider the function evaluation on the right side being performed at the point \( \frac{z_{n+1} + z_n}{2} \). With this, we'll arrive at
\[ B \frac{z_{n+1}-z_n}{\Delta t} = A \frac{z_{n+1} + z_n}{2} \]
The larger down-side of the Crank-Nicolson method is the requirement for solving a large non-symmetric system of equations for every time step.
Leapfrog integration
Leapfrog integration is an explicit symplectic method particularly well suited for mechanical problems, where the accelerations are due to position only, i.e. of the form
\[ \ddot{u}(t) = f(u(t)) \]
which our problem is an example of, as for our problem we have
This might look like Euler's method, but it is subtly different. The position update at \(n+1\) uses the velocity at \(n+1\) and not at \(n\) as for Euler's method.
This still requires solving a system of equations, but now the coefficient matrix is half the size as for the Crank-Nicolson method, and is also symmetric.
Turns out, however, that the mass matrix \(M\) is strongly diagonally dominant and spectrally equivalent to the identity operator. Thus, for large systems \(M\) is quite accurately approximated by a diagonal matrix that has the row (or column) sums of \(M\) as the diagonal elements. The system of equations thus becomes trivial.
