An Arduino-based GPS Disciplined Oscillator

Notes:  August 1, 2016:  The schematic showing the serial interface has been updated.  Its Revision has been changed from Rev. X to Rev. X2.  The circuitry was modified to allow me to use a USB-Serial Adapter Cable manufactured by FTDIChip in lieu of the Elecraft KXUSB cable I had been using.  Please see text below for additional details...

February 6, 2019:  I modified this Frequency-Locked Loop design to incorporate a Phased-Locked Loop based upon Brooks Shera's algorithm (that appeared in the July 1998 issue of QST).  The modifications are described in this blog post:   https://k6jca.blogspot.com/2019/02/an-arduino-version-of-brooks-sheras.html

An Arduino-based GPS Disciplined Oscillator

While I was working on my Antenna Auto-tuner project (see here: k6jca antenna auto-tuner), I discovered that the frequency reading of my HP 8640B at my "second" lab (Lizzie's place in the Californian gold country) was off by a few KHz.

And so, when I finished my tuner, I decided that my next project would be to build a GPS Disciplined-Oscillator (GPSDO) that would then be the master frequency reference for my equipment at Lizzie's.  (The 8640B requires that its External Reference signal be 5 MHz, nominally > 0.5 V peak-to-peak (5V maximum) into 1000 ohms).

Years previously I had built a GPSDO for my home lab, based around a venerable HP 106B Quartz Oscillator that I found, partially disassembled, at a local surplus store.  After repairing it I added Brooks Shera's GPS Frequency Standard (the original QST article and other info can be found here:  http://www.a-aengineering.com/gps.htm) and this has been my home-lab standard for quite some time.

(Click on image to enlarge)

But I certainly didn't have the room for something of this size at Lizzie's, and so I thought, why not make my own, but much smaller, and based around another Arduino NANO processor?

GPSDO Architecture:

Following Brook Shera's QST article, simpler GPSDO designs have been published.  A good example is the GPSDO designed by Bertrand Zauhar, VE2ZAZ, and published in QEX in 2006 (more info here:  http://ve2zaz.net/GPS_Std/GPS_Std.htm). 

This design is based around a frequency-locked loop.  A PIC processor, running from its internal clock, directly counts the 10 MHz frequency-standard oscillator (connected to its T1CKI input) for 16 seconds, and if this count does not equal the "ideal" count value (which is 160,000,000 for a 16 second count period), the PIC adjusts the Oscillator's frequency by changing a DC level created via a PIC-generated PWM signal (filtered with a 2-pole low-pass filter -- note that Brooks Shera used an 18-bit DAC in lieu of PWM).

(Note that, when counting 160,000,000 clocks,  a short counter (e.g. 16 bits or even 8 bits) can be used, so long as the 10 MHz signal is fairly accurate.  If the 10 MHz signal were exactly 10 MHz, a 16-bit counter would read 0x6800 after 16 seconds of counting.  Deviations from this value reflect how far off the actual 10 MHz signal is from its ideal frequency.  An error of 1 count from this value represents an error in frequency of 0.0625 Hz)

The VE2ZAZ design can be simplified even further.  For example, the processor's own clock could be replaced with the 10 MHz oscillator (rather than using a separate clock for this function), and the processor's internal counter set up to use the processor's clock as its clock source, rather than an external source (as was done by VE2ZAZ).  Here are two examples of this architecture:

Murray Greenman, ZL1BPU:  http://www.qsl.net/zl1bpu/MICRO/SIMPLE/SimpleGPS.htm
Hans Summers, G0UPL:  http://www.hanssummers.com/gpsref

The advantage to this method is that some processors limit the frequency of external clocks that clock internal timers.  For example, Arduino resynches external clocks with its processor clock before applying them as clocks to its internal timers. Clock resynchronization within the ATmega328 limits the counter's upper frequency to Fosc / 2.5 (per section 16.3 of Atmel's ATmega328P datasheet).  Assuming the processor is running at 16 MHz (Fosc), then the upper frequency applied to the counter's input should be no more than about 6.4 MHz.

Clearly my 10 MHz oscillator violates Arduino's external-clock upper frequency limit.  So I would need   to divide the 10 MHz by two to create a 5 MHz clock to drive an Arduino timer.  

Or, following in the footsteps of ZL1BPU and G0UPL, if I removed the Arduino's 16 MHz crystal and instead used the 10 MHz oscillator as the processor's clock, I could get around the 6.5 MHz external clock-frequency limitation.  The counter would be clocked with the processor's clock, and, being 10 MHz, this would speed up lock-time by a factor of 2.

But there's a drawback to this latter technique if I want to do more than just discipline an oscillator.  For example, if I want to include serial communications, then there's a problem -- "standard" baud rates invoked via packaged serial-routines would be wrong because the system clock has changed from 16 MHz to 10 MHz, and they would require compensation to accommodate this frequency reduction.  Also, Arduino-specific timing routines, such as millis() and delay() would no longer be accurate, either.

And who knows -- if I removed the 16 MHz crystal, perhaps I would accidentally lift a pcb trace or cause some other damage to the PCB itself.

So I was a bit leery about changing the Nano's processor clock from 16 MHz to 10 MHz.  I decided not to make this hardware mod and, instead, to continue using Arduino's standard 16 MHz processor clock.  And so the oscillator is divided by 2 before applying it, as 5 MHz, to the external-clock input of one of the Arduino timers.

But to which of Arduino's Timers should the 5 MHz be applied?

The choices were Timers 0 and 1 (Timer 2 only uses the processor clock, if I recall correctly).  Timer0 is 8 bits and Timer1 is 16 bits.  But my decision was dependent, in part, on how I generated DC for the oscillator's EFC (Electronic Frequency Control) voltage.

If I used Timer1 to generate a 16-bit PWM signal, filtered to be the DC EFC voltage, then Timer0 would have to be assigned to measuring 5 MHz (or more specifically, the deviation from an ideal 5 MHz count).

But Timer0 is used for Arduino's timing routines, and I didn't know how my proposed use of Timer0 might affect a software UART that I wanted to implement.

But did I really want to use PWM?  Using PWM to generate the EFC voltage was potentially the least-expensive way to create DC (although this design isn't cost-constrained), although it does require that the PWM waveform be filtered to prevent oscillator modulation by the waveform's AC components.

But how much filtering would be required?

Time for a test -- I found a PWM Library (here and more info here)that would let me create a 16-bit PWM signal with Timer1 (pwm.h).

When I tested the routines in this library I discovered that I needed to set the clock frequency to 123 Hz if I were to get a full 16 bits of resolution.  (Note, though, that I didn't do extensive testing, so take this conclusion with a grain of salt).

To prevent the 123 Hz PWM waveform from modulating the 10 MHz oscillator module, I would need to filter it.

Assuming the PWM waveform was a 123 Hz square-wave with an amplitude of 5 Vpp, the amplitude of the fundamental component at 123 Hz would be 1.27 times this, or 3.18 Vpeak. To ensure that this component did not modulate the oscillator, I wanted its amplitude, at the EFC input to the oscillator, to be less than the DAC's lsb step of 0.000153 volts.

Which means it would need to be attenuated by at least 86 dB.  A 4 or 5-pole Butterworth filter with a cutoff frequency of around 10 Hz should be fine.  And it would have a relatively fast settling time, too.

OK, but...wouldn't a 16-bit DAC be much simpler?  There would be less wiring, and, after all, this isn't a commercial product, so I wasn't concerned with eliminating every unnecessary cent of cost by avoiding a relatively more expensive DAC.

In conclusion, I decided to keep the implementation simple.  The EFC voltage is generated with a 16-bit DAC.  Arduino's Timer1 counts the 5 MHz clock, and Timer0 is used for whatever opaque Arduino timing functions it is normally used for.

Here's the basic block diagram of my GPSDO:

(Click on image to enlarge)

Like the VE2ZAZ design, its architecture is based on a frequency-locked loop, not phase-locked loop.

The VCO control voltage (EFC) is generated by a 16-bit DAC, rather than a filtered-PWM signal.

Note that the GPS Receiver's 1 PPS signal drives two inputs.  The first is Timer1's ICP1 input -- the leading edge of the 1 PPS pulse captures a "snapshot" of Timer1's count.  (Timer1 is always running, clocked by 5 MHz.  It never stops).

The leading-edge of the 1 PPS signal also generates an interrupt that signals the software that a count (captured via ICP1) is ready for examination.  (It's possible that ICP1 might also generate an interrupt.  I didn't investigate, as I wasn't concerned with saving I/O pins.  If it does, then INT0 isn't needed.)

GPS Receiver:

Digikey sells a number of different GPS modules.  I chose the GNSS2, from MikroElectronika, because it uses the Quectel L76 GPS receiver.  I'm familiar with this receiver (having used it previously in a commercial product), and so this is the module I decided to go with.  It receives both GPS and GLONASS satellites.

But after I received my GNSS2 module I soon discovered that it had a problem -- specifically, sometimes the L76 would never come out of Reset, and when this occurred, the voltage on the L76's Reset pin would be sitting at about 1.6V instead of 3.3V.

I wasn't driving the board's Reset line.  Instead, I was simply relying on the pullup on GNSS2's RST# line (R5) to keep the L76 Reset signal (T_RST#) high.

Clearly something was wrong, but what?

While probing the module, I noticed that if I touched my scope-probe to the Reset pin of the L76 receiver while it was running, the voltage on this pin would drop from 3.3V to 1.6V and the receiver would lock up.

Not good! 

Per Quectel's L76 Hardware Design Manual, the Reset line should either be left open (no connection) or connected to VCC (note that it doesn't specify if the connection to VCC should be through a pullup or if this pin should be tied directly to the power rail).  And if Reset were to be driven, the manual states that it should be driven with an open-collector driver:

(Click on image to enlarge)

The GNSS2 design, instead of using Quectel's recommended open-collector Reset circuit, drives the L76 Reset line with a TXB0106 Bidirectional Level Shifter (with Auto-sense).

I've had problems with auto-sensing level shifters in the past.  In my opinion they should be used with drivers having well-defined drive levels, not with pullups or with signals whose source circuitry is poorly or incompletely described (such as Quectel's Reset signal).

So, to fix the problem I cut the trace between the TXB0106 and the Reset pin of the L76 and instead connected the L76 Reset pin to one of the unused header-pins on the GNSS2.  I'd have the Nano drive the Reset with Quectel's recommended circuit, just in case, sometime, I needed to Reset the GPS receiver.

And because the Arduino Nano would need access to the GNSS2's 1 PPS signal, I also jumpered this signal from the L76 to another unused header pin on the GNSS2 board.

And finally, note that I changed the PWR SEL jumper on the GNSS2 board from 3.3V to 5V -- the GNSS2 board will be powered from 5V, and its I/O set to 5V, not 3.3V, logic levels.

The mods are shown below, added by me to the MikroElectronika GNSS2 schematic:

(Click on image to enlarge)

And here's a picture showing these mods:

(Click on image to enlarge)

In the picture below, I've connected the GNSS2 module's Serial TX Out to a software-serial port instantiated in the NANO.  The NANO is decoding the GPS NMEA stream and displaying on the LCD the "time received", the "number of satellites used" (15), and the "date" (formatted as February 8, 2016), all decoded from the NMEA stream.

(Click on image to enlarge)

More info on the NMEA stream can be found here:  http://www.quectel.com/UploadImage/Downlad/L76_Protocol_Specification_V1.2.pdf

(I will mention, too, that I'm using the companion MIKROE-363 Active GPS Antenna, also available through Digikey).


10 MHz VCO:

I had a couple of  HP 10811-60111 Quartz Crystal Oscillators in my junkbox.  These are very nice ovenized 10 MHz oscillator modules with an "Electronic Frequency Control" (EFC) input that allows the Oscillator's frequency to be shifted via a voltage applied to this input.  The EFC voltage range is -5 to +5 volts.

To get an idea of how frequency varies with EFC voltage, I tested my two modules (after first adjusting their frequency to be close to 10 MHz when the EFC voltage was forced to be 0 volts):

(Click on image to enlarge)

As you can see from the curves above, maximum dF/dV occurs when the EFC voltage approaches +5 VDC.  At 10 MHz this slope is about -0.19 Hz/V for the newer oscillator (blue line).  It is about -0.29 Hz/V for the older oscillator. 

(In the end, I went with the older oscillator, as it seemed that the newer module might have been experiencing some frequency glitches).



The Schematics:

Arduino NANO schematic page:

Please note that this schematic page has been changed from Rev. X to Rev. X2 (1 August 2016) to incorporate the FTDIChip USB-serial cable, and then from X2 to X3 (9 October 2018) to fix U2 wiring error in schematic (explained below)

(Click on image to enlarge)

Notes on the Arduino NANO schematic page:

• Transistor Q2 is the Quectel-recommended Reset circuit for the L76 GPS Receiver.
• Q1 converts the GPS Receiver's 3.3V 1 PPS signal to 5 volts (and also inverts it so that, as seen by the NANO, the pulse's leading edge is a falling edge).
• I2C bus pullups are on the LCD Module.
• RC networks R12/C18 and R6/C7 provide a bit of low-pass filtering to the 5 MHz and 10 MHz clock edges (I've had past experience with on-board clocks interfering with GPS reception.  Adding high-frequency roll-off (and burying clock traces between PCB layers, an option not open to me in this design) was how we got around this problem in the past).
• U6 is the 16-bit DAC.  Its output is 0-5V.
• U2 converts the DAC's 0-5V output to a range of -5 to +5 V.  I use a TL072 whose power rails are at +/- 12V so that there are no input common-mode range issues (such as I had experienced here).  Rev X3:  Corrected schematic per these changes:  U2.3 connects to U6.1 not U6.6, and U2.5 connects to U6.6 not U6.1 (note that the board was wired correctly).
• U3 provides a "clean" 5V source, separate from the (potentially noisy) 5V_DIG power, as an analog reference voltage (filtered with a single-pole low-pass filter created with R5 and C15).  (Note that U6.1 has a minimum input impedance of 200K ohms, so it is only minimally loading this reference voltage).
• R7 and C12 provide a final single-pole low-pass filter for the EFC voltage, in an attempt to attenuate high-frequency noise generated by the op-amp (see this reference: SLOA049B, specifically figures 13 and 15).
• Header J4 allows me to route the 5 MHz signal to either Timer0's clock input or Timer1's clock input to give me an option as to which timer would be used to count the 5 MHz clock.  (But, per the discussion above, I've jumpered it to select Timer1's input).
• Nano's D10 pin is the Receive-data pin for a software-implemented 9600 baud UART.  It receives the serial NMEA stream from the L76 GPS Receiver.
• Nano's D11 pin is the Transmit-data pin for the same software UART.  This TX signal is inverted by Q3 before going to the "Ring" terminal of J5, a stereo headphone jack (Rev X2: assuming Plug P1 is attached to Header J7), thus allowing me to use an Elecraft KXUSB USB-serial cable to monitor GPSDO performance at any time with my laptop and with a terminal-emulator program, such as puTTY, without removing the GPSDO case covers.  Rev X2:  I can use a less expensive (and more commonly available) FTDIChip USB-serial cable (e.g. from Digikey, p/n TTL-232R-5V-AJ) in lieu of the Elecraft KXUSB by connecting P1 to J6 instead of to J7 (the latter for the Elecraft KXUSB cable).  Doing so removes the Q3 inverter (necessary for the Elecraft cable).  Note that the KXUSB adapter assumes that the "idle" condition (i.e. "Mark") is low and that the Start Bit goes High ("Space"), while the FTDIChip cable assumes that "idle" is High and the Start Bit goes Low.
• R15 helps to "square-up" the 10 MHz duty cycle.  Without it, the duty cycle is about 63/37 instead of the desired 50/50.

Frequency Distribution Schematic page:

(Click on image to enlarge)

Notes on the Frequency Distribution schematic page:

• At the bottom of the page is a DPDT switch.  I use this switch to change the intensity of the LCD's backlight (bright/dim).  It connects to the 2-pin header (normally shorted) on the back of my LCD module.
• VE2ZAZ used an MC3487 for frequency distribution, and I liked the idea, so I did, too.  I did, however, add low-pass filters to each output to knock down the higher harmonics of each clock.
• Note, even though I filter the output clocks: there will still be higher-order harmonics present in the clock waveforms!  After all, I've only added a single-pole filter.  To minimize reflections (and possible problems introduced by the ringing it produces), I recommend using 50 ohm coax between this GPSDO's outputs and any equipment it is connected to.  And these coax cables should be terminated at their far ends with 50 ohms.
• To illustrate the importance of coax termination, here's a screen-shot of my oscilloscope with two 5 MHz signals from my GPSDO connected to its inputs via 6-foot lengths of RG-58A coax.  You can see the ringing in the top trace, where the coax is unterminated at the scope input.  The bottom trace looks much better -- its coax is terminated in 50 ohms at the scope input.

Terminate that Transmission Line!

Power Supply Schematic page:

(Click on image to enlarge)

Notes on the Power Supply schematic page:

• A standard design.  Note that J4 and J6 allow me to connect the front-panel BNC to either the 10811's 10 MHz output or its EFC input (for convenient monitoring of the EFC voltage with external equipment or to let me force it to a fixed level with a DC power supply).
• The 10811-60111 oscillator module requires the following power:  Oven:  20-30 VDC; turn-on load is 42 ohms minimum.  Steady state power will drop to a typical value of 2.0 watts (at 25 degree C and 20 VDC).  Oscillator:  11.0 to 13.5 VDC, 30 mA typical, 40 mA max.  
• 1N4001 diodes could be used in lieu of the 1N4002 diodes I used (they were in my junkbox).  Just keep in mind that whatever diode you use, its PRV rating should be at least twice the secondary's peak voltage.  E.g. the peak voltage of my transformer's secondary (measured on either side of the grounded center-tap) is about 20 V when the primary voltage is 120 VAC.  Assuming the max AC line voltage could be 132 VAC (120 VAC + 10%), then the maximum peak voltage would be 22 V, and so the diode PRV rating should be at least 44 VDC.

The Build

Time to build the circuit and mount it in a chassis.  But what chassis to use?  Here's a swapmeet find I dug out of my "Box o'boxes":

Removing the covers reveals the amount of space available:

The build begins.  The transformer and other sub-assemblies are mounted on a bottom plate fabricated from a piece of copper-clad PCB material that I cut to size.  And note that I try to keep the transformer as far away from the 10811 oscillator module as possible, to minimize the influence of the transformer's magnetic fields on the oscillator.

Final Results:

Here is the back panel.  

A bit tight (because the transformer blocked access to some of the back-panel's area), but everything fit.  Note the SMA for the GPS antenna and the 3.5mm stereo jack in the center for the serial TX data (actually serial RX from the monitoring PC's perspective -- a USB-to-RS232 adapter cable plugs into this jack, thus allowing the user to monitor GPSDO performance).

The front-panel, before gluing on an overlay.

And the front-panel, finished, with its overlay...

The overlay was created using Microsoft Word and then printed on standard printer paper.  I then cut the paper to a rough size slightly larger than the front-panel and glued it to the panel with spray-adhesive (3M Super 77 Adhesive).  I then cut out the holes and painted it with multiple layers of clear Krylon Matte Finish spray paint.

And some holes were drilled in the two side panels to provide a bit of ventilation (you can see them in the image, above). 


Software:

[17 March 2016:  Please note that the software has changed from the original posting in February  (Software revision 160222).  The latest software (Rev. B_160314a)  is described below.]

Measuring Oscillator Frequency:

This GPSDO does not measure frequency.  Rather, it measures deviation from an ideal frequency of 5 MHz.

The 5 MHz clock drives the clock input of the Arduino's Timer1.  This timer never stops -- it is always being clocked, and thus always incrementing with each clock pulse, until it reaches 0xFFFF, rolls over to 0x0000, and continues climbing again.

Once a second, triggered by the leading edge of the GPS Receiver's 1 PPS signal, the software takes a snapshot of Timer1's count.  It then calculates the delta between this latest snapshot and the snapshot it took a second earlier.

If the counter were being clocked at exactly 5 MHz (and if there were no jitter on the 1 pps signal), this delta should always be 19264 (decimal).

If this delta is less than 19264, then the oscillator frequency is too low.

And if this delta is greater than 19264, then the oscillator frequency is too high.

Every second, software is checking if the delta is too high, too low, or spot-on. 

As we train the oscillator closer and closer to its ideal 10 MHz frequency, the snapshot period must become longer and longer to increase our measurement resolution. That is, if I increase my count time from one second to two seconds, I've increased my counter's "frequency resolution" by a factor of 2.

So, Timer1 has an effective resolution of 1 Hz with a 1-second snapshot period, 0.5 Hz with a 2-second snapshot period, and 0.25 Hz with a 4-second snapshot period.

My original software counted for a fixed duration (in seconds) and then checked, when that duration was finished, if the majority of deltas had been high, low, or unchanged.

Duration-intervals were adaptive in power of two -- the minimum interval was 16 seconds, and the maximum interval was 32768 seconds, providing counter-resolutions from 0.0626 Hz to 30.5 microHertz, per the table below.

(Click on image to enlarge)

But there's a potential problem with this method -- if, for some reason, the oscillator frequency begins to deviate during a long interval, the software won't catch that deviation until the end of the interval, which means that a significant frequency error could occur.

I tried some methods to jump out of an interval if the error grew too large while it was "integrating", but I felt like I was putting band-aids on the algorithm, and so I decided to take another tack...

I noticed that as frequency began to deviate from the ideal, I might see a delta value of 19263 followed by 19265 (with perhaps some 19264 values in between), or the opposite: 19265 followed by 19263.  These low/high (or high/low) pairs would occur more frequently as the frequency deviated, and I wondered what would happen if, instead of having a fixed count-duration (as described above), I instead let the counter count until I saw two 19263 or 19265 values in a row (19264 values between the two would be allowed).

It turns out, that although this new algorithm is not perfect, it works pretty well. 

The New Algorithm:

There are no longer intervals of fixed durations. Instead, the counting-interval is variable, and it is based on the following procedure:

1.  Initialize a "seconds" counter to 0.

2.  Once a second, increment the seconds-counter and:

     o  Take a snapshot of the 5 MHz counter and calculate the
         delta from the snapshot of the previous second.

     o  If this delta and the last non-19264 delta are both 
         greater than19264, the frequency has changed.  
         Record the elapsed seconds and go to step 3.

     o  If this delta and the last non-19264 delta are both 
         less than19264, the frequency has changed.  
         Record the elapsed seconds and go to step 3.

      o  Otherwise, repeat step 2.

3.  Using the elapsed-seconds reading from step 2, calculate the frequency error and use this value to determine the amount the DAC should be changed to correct for this error.  Update the DAC.

4.  Return to Step 1.

Note that the calculation to determine the DAC correction is:

(1 cycle(error) / elapsed_seconds) / ((Osc_Hz/V @ 5 MHz) * (DAC_V/step))

From the table, above:

• Osc_Hz/V @ 5 MHz = 0.145 Hz/V.
• DAC_V/step = 1.53E-4 V/step

Plugging these values into the equation above, the DAC correction for this oscillator is:

DAC_Correction = 45075 / elapsed_seconds

Note, my code actually uses 42240 instead of 45075, as I wanted to ensure that the overall "gain" was less than 1 (42240 = 15 * 2816, using the 2816 calculated from my earlier algorithm for a 16-second interval and scaling this to "almost" 1 second by using 15 in lieu of 16).

This algorithm works surprisingly well for short runs of elapsed-seconds.  But I found for longer periods, oscillator drift (i.e. aging) cannot be compensated for simply with the algorithm above.

Below is a plot of frequency-error, and its correction over time.  Please note:

• This plot shows the frequency deviation between my Arduino-based frequency-locked GPSDO and my older HP 106B-based phase-locked GPSDO (the latter used the frequency "reference" for these measurements.
• An HP 5370B Universal Time Interval Counter measured the time-interval between the rising edge of the Arduino-based GPSDO's 5 MHz output (applied to the counter's START input) and the rising edge of the 106B-based GPSDO's 5 MHz output (applied to the STOP input).
• In the plot below, this time-interval (the blue line) was measured once every 60 seconds.  The amber line is the frequency error, in parts-per-billion.  It is the slope of the blue line, taken between two consecutive 60-second time-interval samples.

(Click on image to enlarge)

(Matlab tool written by Dick Benson, W1QG)

Even though this plot was generated with my HP 106B-based GPSDO as a frequency reference (which will have its own frequency errors -- a Cesium standard would have been much nicer!), the error between the two oscillators is typically between 0.1 ppb and 0.01 ppb.

That is, worse-case error is about 1 Hz at 10 GHz.  Not too bad!

(An interesting note -- per the blue-line in the plot, with this algorithm the clock-difference delta between DAC updates is typically about 200 ns (i.e. one 5 MHz clock cycle).  Sometimes this difference is 400 ns, though.  I'm not sure why it isn't always 200 ns, but...it isn't.)

Note that clock-difference (blue line) is always increasing.  This is because the frequency of my 10811 oscillator module is always increasing (note that HP specs its Aging Rate at < 5E-10 / day).

I wondered -- could I compensate for this aging and improve accuracy?

Ideally, such an algorithm would be an adaptive one, automatically determining the aging-rate and adapting the GPSDO's compensation as the rate changes over time.

But I couldn't find a good way to do this, so I took a simpler, more brute-force approach:  simply measure the aging rate externally (with my HP 5370B) and plug this value into the code.

To measure the aging rate I fixed the DAC value at a value that gave me, essentially, almost exactly 5 MHz, and then I measured how this frequency changed over time.  Here's that plot:

(Click on image to enlarge)

The error deviation (slope of the amber line) is about 0.25 ppb in 12 hours, or 0.02 ppb per hour (i.e. right at HP's Aging Rate spec of < 0.5 ppb per day).

So, assuming a fixed Aging Rate of 0.02 ppb per hour and knowing that my DAC correction-factor is 0.0044 ppb per DAC step, I should apply a correction of about 4.5 DAC steps per hour of "elapsed-time" that the algorithm measured while looking for a frequency change.

Here's plot with Aging-Rate compensation included in the algorithm:

(Click on image to enlarge)

If you compare this plot to the previous plot, you'll see that this new algorithm's frequency-error is below 0.01 ppb more often. BUT, notice the glitch at about 30 hours -- error was almost 0.2 ppb.  This spike in the frequency-error occured because, for some reason, the slope of the blue-line (the clock-difference plot) was flat over the 5-hour period from about 20 hours to 25 hours -- that is, there wasn't any noticable aging.  But the algorithm didn't know this -- it has no way to determine if aging is occuring, or not.  It simply assumes that aging is always occuring, and thus it assumed that the oscillator had been aging for those 5 hours, too.

Therefore, when the algorithm determined that the 5 MHz frequency had finally changed (after about 10 hours (!!!) of elapsed-time, around the 30-hour mark), it overcompensated, which resulted in a frequency change that was larger than it should have been.

Why was the clock-difference flat for about 5 hours?  I don't know.  But it points to a potential flaw in this method of compensation.  If the Aging Rate is not constant, it's possible that there will be brief periods where the frequency-error will exceed 0.1 ppb as it compensates for the assumed constant rate.

So is this compensation worth pursuing? Frankly, for my application, not really, given that error still drifts up towards 0.1 ppb and that now there's the possibility of an over-correction.  But perhaps the algorithm could be improved upon, which is why I include the discussion here.

(A note on the plots above:  I believe the results above are fairly representative of the accuracy of the Arduino-based GPSDO.  The 106B-based GPSDO is a phase-locked architecture and thus much different from the frequency-locked architecture of the Arduino-based GPSDO.  And I believe these architecture differences are significant enough to ensure that one GPSDO isn't simply mimicing the other GPSDO and thus giving an erroneous picture of frequency accuracy).


Updating the DAC:

When within the 1-second window between 1 PPS pulses should the DAC be updated and the LCD written-to?

One might think, just do these tasks after the 1 PPS pulse occurs, when you've taken the snapshot of the count.

But there's a potential problem.  I'll explain using the image below:

The top trace is the NMEA stream and the bottom trace is the 1 PPS pulse from the L76.

As you can see, in the image above, the 1 PPS pulse occurs while data is being received.

I noticed that the software UART could miss NMEA data if Arduino code was, for example, writing to the LCD at the same time the NMEA stream was being received.

So, clearly, the best time to write to the LCD would be after the NMEA packet has finished.

No problem -- I can easily detect the end of a packet.

But, if a packet is short, then the 1 PPS pulse can occur after the end of a packet (rather than in the middle of a packet, as shown above).

So the method is this:  the delta is calculated, the votes are summed, decisions are made, and the LCD is written to, if the NMEA packet has ended and if I've seen the leading edge of the 1 PPS pulse.

This results in an interesting display artifact resulting from the variation in time between the end of the NMEA packet and the leading-edge of the 1 PPS pulse.  Specifically, the "seconds" field in the time display on my LCD doesn't always update "on the beat."  Sometimes it is slightly off.  Not a big deal for me, but worth mentioning).

After the DAC is updated I include a one-second "settling time" in my state-machine to ensure that the DAC has stabilized before I start counting seconds and accumulating votes.


Software Serial Port:

I wanted to be able to examine the serial-data being streamed by the GPS receiver and extract from it useful information such as date, time, and the number of satellites used in the position-fix.

I also wanted to be able to send, from the GPSDO, a serial-data stream (in CSV format) containing GPSDO performance data that could be captured via my laptop and easily analyzed using, say EXCEL.

The Arduino has a hardware serial port, but this port is used for downloading code or uploading, for example, "Serial.print" data generated via the software running on the Arduino.  But I wanted to simultaneously receive the NMEA stream from the GPS receiver and also transmit, via an external serial port, performance data, even if the Arduino's hardware serial port were in use with other tasks (e.g. debug tasks).

So I used the SoftwareSerial library to implement a software serial port (instantiated with the name mySerial) running at 9600 baud (the GPS Receiver's default baud rate).  The RX side of this port receives the NMEA data from the GPS receiver, while the TX side sends "performance data" out via a connector on the GPSDO back panel to my laptop.


NMEA Decoding:

To keep NMEA decoding from delaying other software processes (which could happen if I first received the entire NMEA packet before attempting to decode it), I decode the relatively-slow NMEA stream on a character-by-character basis ("on the fly"), as each new character is received, using a simple state-machine that tracks where characters are in the NMEA stream and then copies, as data arrives, the appropriate information (time, date, satellites-used) for later display on the LCD.


LCD Display:

[With the latest revision of the code the data displayed on the LCD has changed slightly. The description below is for the previous revision of the software, which is close.] 

In addition to supplying useful information such as time and date, the LCD display also gives me GPSDO status at-a-glance.

The fields in the LCD display above are:

Top Line:

• "15:59:11" -- Time (hh:mm:ss format) extracted from the GPS NMEA serial-data stream.
• "SV: 18"  -- Number of satellites used for the Position Fix (18 satellites in this instance).  Also from the NMEA stream.  Indicates how well the GPS receiver and its antenna are working.

Bottom Line:

• "02/21/16" -- Date string (converted to mm/dd/yy format) from GPS NMEA stream.
• "83DE" -- DAC value, in Hex (i.e. 0x83DE.  Total range is 0x0000 - 0xFFFF).
• "3" -- This value represents the GPSDO's current Integration Level (4096 second of integration for Level 3), .  The Integration Level is represented by a single digit, thus, the two top integration levels (10 and 11 at this time) are represented by characters A and B.  Lower levels use numbers, such as "3" for Level 3, as seen in the display above.  [This character is no longer displayed with the latest revision of the code.]
• "+" -- The "+" signifies that the "votes" accumulated during the previous integration-period resulted in an increase in the DAC value.  A decrease is represented by a "-", and no-change is represented by "."

The Code:

Below is a complete listing of my Arduino code.  Take it with a block of salt, as I consider myself much more of a hardware designer than a software coder (and if you have any recommendations or corrections, please feel free to write me!).

Also, this blog website doesn't seem to have a simple way to add code to posts, so I've cobbled up the box below.  Note that it will wrap lines around, which makes reading difficult.  So my recommendation is that you just do a simple copy of the code in the box below and paste it into a text editor that you prefer.

//------------------------------------------------
// Author:  Jeff Anderson, K6JCA
//
// Rev:  B_160314a
//
//------------------------------------------------
#include  <LiquidCrystal_I2C.h> 
#include  <SoftwareSerial.h> 
#include  <PWM.h> 
#include  <Wire.h> 
#include  <stdlib.h> 


// Compiler Directives
//#define ECHO_NMEA
//#define PRINT_GPSDO_STATES
//#define PRINT_ACCUM_TOTAL
#define PRINT_CSV
//#define PRINT_DELTA
//#define HOLD_DAC

// Constant Defines
#define MAX5217  0x1C  // DAC I2C address
#define INT0_PIN 2     // NANO pin, INTERRUPT 0
#define GPS_RST  7     // NANO pin, reset GNSS-2 receiver
#define TX_IN   10     // NANO pin, software UART RX
#define ICP1     8     // NANO pin, T1's ICP1 input
#define T1_CLK   5     // NANO pin, T1's clock input

// Starting Values
#define MAX_DAC_STEP     2816   // 16-second DAC step
#define GAIN             15.0
//#define GAIN             0.0
#define START_DAC        0x8000

#define OSC_PPB_DRIFT_PER_HOUR      0.02   // measured osc aging 
//#define OSC_PPB_DRIFT_PER_HOUR      0.0     // DON'T CORRECT FOR AGING
#define PPB_CORRECTION_PER_DAC_STEP 0.0044  // from osc. measurements

#define IDEAL_DELTA 19264   // for a 1 second period.

// gpsdo states:
#define WAIT   0
#define START  1
#define SETTLE 2
#define GO     3
#define SETTLE_2 4

// set the LCD address to 0x27
LiquidCrystal_I2C lcd(0x27,16,2); 
// Define software UART, Rx pin 10, Tx pin 11
SoftwareSerial mySerial(10, 11);  


int byteGPS;
char gps_time[9];
char gps_date[9];
char sats_used[3];      

int nmea_decode_state;   // state of NMEA decode process (to extract date, time, and SV)

boolean date_flag;
boolean time_flag;
boolean sv_flag;
boolean lcd_update_flag;
boolean first_lcd_update;
boolean nmea_start_flag;  // signals start of nmea packet
boolean nmea_end_flag;    // signals end of nmea packet
boolean one_pps_flag;     // signals leading-edge of one_pps_pulse

unsigned int prev_t1_count;
unsigned int t1_count;
unsigned int t1_delta;

int number_of_sats;

byte gpsdo_state;

unsigned int dac;
long         dac_delta;
unsigned long accum_seconds;
unsigned long long_dac;

int  dac_direction;   // +1 = positive change, -1 = neg. change, 0 = no change.

char tick_count;
boolean display_countdown;
unsigned int prev_countdown;

unsigned long zero_count;
unsigned long prev_count;
boolean one_flag;
boolean minus_one_flag;

float float_zero_count;

int  last_direction;
long run_error;
long accum_error;
long prev_accum_error;

float drift_corr_per_hour;  // drift correction, per hour

void setup()
{
  
  Serial.begin(57600);  // init serial port
  lcd.init();           // initialize the lcd 
  InitTimersSafe();     //initialize all timers except for 0, to save time keeping functions

  pinMode(INT0_PIN, INPUT);
  pinMode(TX_IN, INPUT);
  pinMode(ICP1, INPUT);
  pinMode(T1_CLK, INPUT);
  pinMode(GPS_RST, OUTPUT);

  // Reset for 15 msec
  digitalWrite(GPS_RST,HIGH); 
  delay(15);
  digitalWrite(GPS_RST,LOW); 
  
  date_flag = false;
  time_flag = false;
  sv_flag = false;
  lcd_update_flag = false;

  nmea_start_flag = false;
  nmea_end_flag = false;
  one_pps_flag = false;

  first_lcd_update = true;

  nmea_decode_state = 0;

  number_of_sats = 0;

  tick_count = 0;
  display_countdown = false;
  prev_countdown = 0;


  zero_count = 0;
  accum_seconds = 0;
  dac                 = START_DAC;
  dac_direction       = 0;
  gpsdo_state = WAIT;

  dac_delta = 0;
  long_dac = long(dac);

  one_flag = false;
  minus_one_flag = false;

  
  last_direction = 0;
  run_error = 0;
  prev_count = 0;
  accum_error = 0;

  drift_corr_per_hour = OSC_PPB_DRIFT_PER_HOUR/PPB_CORRECTION_PER_DAC_STEP;
  
  write_to_dac(dac);

  //  Print a message to the LCD...
  //
  lcd.backlight();
  lcd.print("  Hi, I'm your");
  lcd.setCursor(0,1);
  lcd.print("  K6JCA GPSDO!");

  mySerial.begin(9600);
  mySerial.println("Hello Jeff");
  
  delay(2000);
  lcd.clear();
  lcd.print("  Waiting for   ");
  lcd.setCursor(0,1);
  lcd.print(" NMEA stream... "); 


  // interrupt on falling edge of one_pps signal.
  attachInterrupt(digitalPinToInterrupt(INT0_PIN), one_pps, FALLING); 

  // Timer 1: external clock.  Neg. Edge capture.  No noise cancellation.
  timer1_setup (0x00, -1, 0x00, 0x00, 0x00); 


  #ifdef PRINT_CSV
    Serial.println("Total Seconds,DAC,DacDirection,ZeroCount,RunError,PrevAccumError,DacDelta,New-AccumError,SV");
  #endif 
    mySerial.println("Total Seconds,DAC,DacDirection,ZeroCount,RunError,PrevAccumError,DacDelta,New-AccumError,SV");

}

void loop() 
{
  // Read a byte of the serial port, 
  // (returns -1 if no byte, I believe).
  // Byte available.  Get it and decode NMEA time,
  // date, and SV info.
  byteGPS = mySerial.read();        
  if (byteGPS != -1) {            
    #ifdef ECHO_NMEA
      Serial.write(byteGPS);
    #endif 
    decode_nmea();
  }
  // only update LCD (and do other calcs) 
  // if not receiving NMEA packet AND 
  // if 1pps pulse has happened.
  if (nmea_end_flag && one_pps_flag) {
    lcd_update_flag = true;
    nmea_start_flag = false;
    nmea_end_flag = false;
    one_pps_flag = false;
  }
  if (lcd_update_flag) {
    // convert ascii to number
    number_of_sats = atoi(sats_used); 
    if (number_of_sats >= 4) {
      // 1 pps seems to be present if satellites >= 4 
      // so call the gpsdo state machine.    
      gpsdo_state_machine(); 
    }
    update_lcd();
    lcd_update_flag = false;
    #ifdef PRINT_GPSDO_STATES
      Serial.println(t1_delta);
      Serial.print("state: ");
      Serial.print(gpsdo_state);
    #endif
  }
  
}


void decode_nmea()
{
  // The NMEA packet from the Quectel L76 receiver starts
  // with $GNRMC and ends with $GNGLL, and
  // they arrive in the following order:
  //    $GNRMC
  //    $GPVTG
  //    $GPGGA
  //    $GNGSA (might be more than one GNGSA line)
  //    $GPGSV (might be more than one line)
  //    $GLGSC (might be more than one line)
  //    $GNGLL
  //
  // From this packet of information, I want to:
  //   o  Decode 3 numbers from the gps receiver's
  //      NMEA stream:
  //       o  Time
  //       o  Date
  //       o  SV (number of satellites used for fix)
  //   o  Flag when the packet has started, and also 
  //      flag when it ends.
  //
  //  Extract Time and Date from the $GxRMC command 
  //  (decode states start at state 20).
  //  This command will also set the nmea_start_flag.
  //
  //  Extract SV from the $GnGGA command 
  //  (decode states start at state 3)
  //
  //  When we have all the data we want, we can write
  //  to the LCD, but don't do this until the entire 
  //  NMEA packet has been received. 
  //  Otherwise, writing to the lcd
  //  might cause NEMA info to be missed 
  //  (and not echoed to serial port).
  //
  //  The last line sent in a Quectel L76 NMEA 
  //  packet seems to be the $GNGLL line 
  //  (the first is the $GxRMC line. 
  //  The line-feed is the last character in that 
  //  line, and therefore the last character in 
  //  the data packet. 
  //  (Decode of packet end starts at state 50).
  //  Otherwise, return false if haven't yet 
  //reached the end of packet. 

  switch (nmea_decode_state) {
    case 0:
      if (byteGPS == 'G') nmea_decode_state = 1;
      else {
        if (byteGPS == 'R') nmea_decode_state = 20;
        else nmea_decode_state = 0;
      }
    break;

    case 1:
      if (byteGPS == 'G') nmea_decode_state = 2;
      else {
        if (byteGPS == 'L') nmea_decode_state = 50;
        else nmea_decode_state = 0;
      }   
    break;
   
    case 2:
      if (byteGPS == 'A') nmea_decode_state = 3;
      else nmea_decode_state = 0;
    break;

    case 3:   // 'GGA' found, 
      // now start counting commas to the SV field
      if (byteGPS == ',') nmea_decode_state = 4;
      else nmea_decode_state = 0;
    break;

    case 4:   // found one comma
      if (byteGPS == ',') nmea_decode_state = 5;
    break;

    case 5:   // found 2 commas
      if (byteGPS == ',') nmea_decode_state = 6;
    break;

    case 6:   // found 3 commas
      if (byteGPS == ',') nmea_decode_state = 7;
    break;

    case 7:   // found 4 commas
      if (byteGPS == ',') nmea_decode_state = 8;
    break;

    case 8:   // found 5 commas
      if (byteGPS == ',') nmea_decode_state = 9;
    break;

    case 9:   // found 6 commas, look for 7th...
      if (byteGPS == ',') nmea_decode_state = 10;
    break;

    case 10: 
      // found 7 commas, 
      // next characters should be SV field...
      sats_used[0] = byteGPS;
      nmea_decode_state = 11;
    break;

    case 11:
      if (byteGPS != ',') {
        sats_used[1] = byteGPS;
        sats_used[2] = '\0';
      }
      else sats_used[1] = '\0';
      sv_flag = true;
      nmea_decode_state = 0;
    break;

    case 20: 
      // found an 'R'. 
      // Now let's check if the next characters 
      //are M and C...
      if (byteGPS == 'M') nmea_decode_state = 21;
      else nmea_decode_state = 0;
    break;

    case 21: 
      // found an 'M'. 
      //Now let's check if the next character is C...
      if (byteGPS == 'C') {
        nmea_decode_state = 22;
        // found 'RMC' -- so set nmea_start_flag
        nmea_start_flag = true;
      }
      else nmea_decode_state = 0;
    break;

    case 22: 
      // found a 'C'. 
      // Now let's check if the next character 
      // is a comma
      if (byteGPS == ',') nmea_decode_state = 23;
      else nmea_decode_state = 0;
    break;

    case 23:  // found a comma. 
      //Next characters will be time...
      gps_time[0] = byteGPS;
      nmea_decode_state = 24;
    break;

    case 24:  // time...
      gps_time[1] = byteGPS;
      gps_time[2] = ':';
      nmea_decode_state = 25;
    break;

    case 25:  // time...
      gps_time[3] = byteGPS;
      nmea_decode_state = 26;
    break;

    case 26:  // time...
      gps_time[4] = byteGPS;
      gps_time[5] = ':';
      nmea_decode_state = 27;
    break;

    case 27:  // time...
      gps_time[6] = byteGPS;
      nmea_decode_state = 28;
    break;

    case 28:  // time...
      gps_time[7] = byteGPS;
      gps_time[8] = '\0';
      nmea_decode_state = 29;
      time_flag = true;
    break;

    case 29: 
      // Now let's start counting commas to 
      // the date field (there should be 8)
      if (byteGPS == ',') nmea_decode_state = 30;
    break;

    case 30: 
      // found a 1st comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 31;
    break;

    case 31: 
      // found a 2nd comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 32;
    break;

    case 32: 
      // found a 3rd comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 33;
    break;

    case 33: 
      // found a 4th comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 34;
    break;

    case 34: 
      // found a 5th comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 35;
    break;

    case 35: 
      // found a 6th comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 36;
    break;

    case 36: 
      // found a 7th comma. 
      // check if the next character is a comma
      if (byteGPS == ',') nmea_decode_state = 37;
    break;

    case 37: 
      // found 8 commas. 
      // Now get the date. It's in ddmmyy format, 
      // but I'll display it as mm/dd/yy
      gps_date[3] = byteGPS;
      nmea_decode_state = 38;
    break;

    case 38: 
      // found 8 commas.  Now get the date
      gps_date[4] = byteGPS;
      gps_date[5] = '/';
      nmea_decode_state = 39;
    break;

    case 39: 
      // found 8 commas.  Now get the date
      gps_date[0] = byteGPS;
      nmea_decode_state = 40;
    break;

    case 40: 
      // found 8 commas.  Now get the date
      gps_date[1] = byteGPS;
      gps_date[2] = '/';
      nmea_decode_state = 41;
    break;

    case 41: 
      // found 8 commas.  Now get the date
      gps_date[6] = byteGPS;
      nmea_decode_state = 42;
    break;

    case 42: 
      // found 8 commas.  Now get the date
      gps_date[7] = byteGPS;
      gps_date[8] = '\0';
      date_flag = true;
      nmea_decode_state = 0;
    break;


    case 50: 
      // 'GL' found, now see if there's 
      // another 'L'...
      if (byteGPS == 'L') {
        nmea_decode_state = 51;
        //Serial.println("s");
      }
      else nmea_decode_state = 0;
    break;

    case 51: 
      // 'GLL' found. 
      // Now look for end of $GnGLL 
      // packet (Line Feed).
      if (byteGPS == 10)  {
        nmea_end_flag = true;   
        nmea_decode_state = 0;   
      }
    break;
      
    default:
      nmea_decode_state = 0;
    break;
   
  }
}

void update_lcd()
{
  // first, get magnitude of run_error and
  // magnitude of accum_error
  long mag_run_error;
  long mag_accum_error;
  if (run_error < 0) mag_run_error = -run_error;
  else mag_run_error = run_error;
  if (accum_error < 0) mag_accum_error = -accum_error;
  else mag_accum_error = accum_error;
  
  // first time into this call, clear the screen. 
  // Otherwise, just write w/o clearing.
  if (first_lcd_update == true) {
    lcd.clear();
    first_lcd_update = false;
  }
  tick_count++;
  if (tick_count == 6) { 
    tick_count = 0;
  }

  lcd.setCursor(0,0);
  lcd.print(gps_time);
  lcd.print("  ");
  lcd.setCursor(10,0);
  lcd.print("SV: ");
  // compensate if number of sats < 10.
  if(sats_used[1] == '\0') lcd.print(" ");
  lcd.print(sats_used);
  lcd.setCursor(0,1);
  lcd.print(gps_date); 
  lcd.print(" ");
  lcd.setCursor(9,1);

  switch(tick_count) {
    case 2:
      // Print current value Zero count
      if (zero_count < 10000) lcd.print(" ");
      if (zero_count < 1000) lcd.print(" ");
      if (zero_count < 100) lcd.print(" ");
      if (zero_count < 10) lcd.print(" ");
      lcd.print(zero_count);
    break;

    case 3:
      // Print previous zero-count
      if (prev_count < 10000) lcd.print(" ");
      if (prev_count < 1000) lcd.print(" ");
      if (prev_count < 100) lcd.print(" ");
      if (prev_count < 10) lcd.print(" ");
      lcd.print(prev_count);
    break;

    case 4:
      // Print magnitude of previous run-error.
      lcd.print("rER");
      if (mag_run_error < 1000) lcd.print(" ");
      if (mag_run_error < 100) lcd.print(" ");
      if (mag_run_error < 10) lcd.print(" ");
      lcd.print(mag_run_error);
    break;

    case 5:
      // Print magnitude of current accumulated_error.
      lcd.print("aER");
      if (mag_accum_error < 1000) lcd.print(" ");
      if (mag_accum_error < 100) lcd.print(" ");
      if (mag_accum_error < 10) lcd.print(" ");
      lcd.print(mag_accum_error);
    break;

    default:
      // Print DAC value and sign of last delta.
      // print leading zeroes (normally suppressed
      lcd.print(" ");
      if (dac < 4096) lcd.print("0");
      if (dac < 256) lcd.print("0");
      if (dac < 16) lcd.print("0");
      lcd.print(dac,HEX);
      lcd.print(" ");
      dac_direction_to_lcd();
  
    break;
  }
}

void one_pps()
{
  one_pps_flag = true;
}

void timer1_setup (byte mode, int prescale, byte outmode_A, byte outmode_B, byte capture_mode)
{
  // NOTE:  This code found at:
  // http://sphinx.mythic-beasts.com/~markt/ATmega-timers.html
  
  // enforce field widths for sanity
  mode &= 15 ;
  outmode_A &= 3 ;
  outmode_B &= 3 ;
  capture_mode &= 3 ;

  byte clock_mode = 0 ; // 0 means no clocking - the counter is frozen.
  switch (prescale)
  {
    case 1: clock_mode = 1 ; break ;
    case 8: clock_mode = 2 ; break ;
    case 64: clock_mode = 3 ; break ;
    case 256: clock_mode = 4 ; break ;
    case 1024: clock_mode = 5 ; break ;
    default:
      if (prescale < 0)
        clock_mode = 7 ; // external clock
  }
  TCCR1A = (outmode_A << 6) | (outmode_B << 4) | (mode & 3) ;
  TCCR1B = (capture_mode << 6) | ((mode & 0xC) << 1) | clock_mode ;
}

void gpsdo_state_machine()
{
  switch (gpsdo_state) {
   
    case WAIT:
      gpsdo_state = START;
    break;

    case START:
      #ifdef PRINT_CSV
        Serial.print(accum_seconds);
        Serial.print(",");
        Serial.print(dac);
        Serial.print(",");
      #endif 
      mySerial.print(accum_seconds);
      mySerial.print(",");
      mySerial.print(dac);
      mySerial.print(",");
     
      write_to_dac(dac);
      gpsdo_state = SETTLE;
      accum_seconds++;
    break;

    case SETTLE:
      // DAC has been set.  Let's let it settle.
      gpsdo_state = SETTLE_2;
      accum_seconds++;
    break;

    case SETTLE_2:
      // DAC has been set.  Let's let it settle.
      // Meanwhile, get current count and 
      // set to previous count.
      t1_count = ICR1;
      prev_t1_count = t1_count; 
      gpsdo_state = GO;
      accum_seconds++;
    break;

    case GO:
      accum_seconds++;  // another second has passed
      prev_t1_count = t1_count; // prev. counter snapshot
      t1_count = ICR1;          // new counter snapshot

      // compensate for counter wraparound
      if (t1_count < prev_t1_count) {
        t1_delta = 65536 - (prev_t1_count - t1_count);
      }
      else {
        t1_delta = t1_count - prev_t1_count;
      }
      #ifdef PRINT_DELTA
        Serial.print(t1_count);
        Serial.print(" - ");
        Serial.print(prev_t1_count);
        Serial.print(" = ");
        Serial.println(t1_delta);
      #endif
      // a positive delta following a negative delta
      // cancel. So run of zeroes is not yet over.
      // continue counting zeroes!
      if (t1_delta > IDEAL_DELTA) {
        if (minus_one_flag) {
          // +1 and -1 cancel so
          // continue counting
          dac_direction = 0;
          zero_count = zero_count + 1;
          minus_one_flag = false;
        }
        else {
          if (!one_flag) {
            // Previous delta was 0, set set flag
            // that this delta is a 1.
            // Don't stop counting zeroes yet.
            one_flag = true;
            zero_count = zero_count + 1;
            dac_direction = 0; 
          }
          else {
            // Two positive deltas in a row.
            // Run of zeroes is over!
            // Reset flag and indicate direction to 
            // change dac.
            // Positive deltas mean increase dac
            // to lower the oscillator frequency.
            dac_direction = +1;
            one_flag = false;
          }
        }
      }
      else {
        // same idea as above, but instead with a 
        // negative deltampreceeding a pos. delta.
        if (t1_delta < IDEAL_DELTA) {
          if (one_flag) {
            // +1 and -1 cancel, so continue with zeroes
            dac_direction = 0;
            zero_count = zero_count + 1;
            one_flag = false;
          }
          else {
            if (!minus_one_flag) {
              minus_one_flag = true;
              zero_count = zero_count + 1;
              dac_direction = 0;
            }
            else {
              // two minus deltas in a row.
              // Run of zeroes is over!  Reset flag
              // and identify direction to move DAC.
              // Negative deltas mean decrease dac
              // to raise the oscillator frequency.
              dac_direction = -1;
              minus_one_flag = false;
            }
          }
        }
        else {  // the run of zeroes continues...
          zero_count++;
        }
      }
      if (dac_direction != 0) {
        // Run of zeroes has ended.  From run length
        // determine dac correction factor (= run_error).
        // Will also add to that value a "drift" offset
        // (proportional to run length)
        prev_count = zero_count;  // for lcd display purposes
        float_zero_count = float(zero_count);
        prev_accum_error = accum_error;  // for print purposes
        accum_error = long(drift_corr_per_hour * float_zero_count / 3600);
        run_error = long(MAX_DAC_STEP * (GAIN / (float_zero_count + 1)));  // add one so don't divide by 0.
        if (dac_direction == -1) {
          run_error = -run_error;
          last_direction = -1;  // For lcd display purposes
        }
        else {
          last_direction = +1;
        }
        dac_delta = run_error + accum_error;

        #ifdef PRINT_CSV
          Serial.print(dac_direction);
          Serial.print(",");
          Serial.print(zero_count);
          Serial.print(",");
          Serial.print(run_error);
          Serial.print(",");
          Serial.print(prev_accum_error);
          Serial.print(",");
          Serial.print(dac_delta);
          Serial.print(",");
          Serial.print(accum_error);
          Serial.print(",");
          Serial.println(number_of_sats);
        #endif
        mySerial.print(dac_direction);
        mySerial.print(",");
        mySerial.print(zero_count);
        mySerial.print(",");
        mySerial.print(run_error);
        mySerial.print(",");
        mySerial.print(prev_accum_error);
        mySerial.print(",");
        mySerial.print(dac_delta);
        mySerial.print(",");
        mySerial.print(accum_error);
        mySerial.print(",");
        mySerial.println(number_of_sats);

        // Calculate new dac value.
       
        #ifdef HOLD_DAC
          dac_delta = 0;
        #endif
       
        long_dac = long(dac) + dac_delta;
        if (long_dac > 65535) long_dac = 65535;
        if (long_dac < 0) long_dac = 0;
        dac = uint16_t(long_dac);

        // dac updated.  Start counting 0's anew...
        zero_count = 0;
        dac_direction = 0;
        gpsdo_state = START;
      }
    break;

    default:
    break;
  }
}

void write_to_dac(unsigned int dac_value)
{
  // writing 16-bit word to Maxim MAX5217
  
  byte msbyte;
  byte lsbyte;

  msbyte = byte((dac_value & 0xFF00)>>8);
  lsbyte = byte(dac_value & 0x00FF);
  Wire.beginTransmission(MAX5217);
  Wire.write(0x01);  // control word
  Wire.write(msbyte);  // MS Byte
  Wire.write(lsbyte);  // LS Byte
  Wire.endTransmission();
}


void dac_direction_to_lcd()
{
  if (last_direction == +1) lcd.print("+");
  else {
    if (last_direction == -1) lcd.print("-");
    else lcd.print(".");
  }
}

Please note that this code includes Aging-Rate compensation (OSC_PPB_DRIFT_PER_HOUR).


Additional Testing:

EFC noise:

Noise on the EFC will translate into phase noise (or FM modulation) on the 10 MHz signal.  So it's important that the EFC noise be minimized.

To measure EFC noise, I connected the EFC voltage (accessed via my GPSDO's front-panel BNC) to a SigLab 50-21.  Looking at the spectrum from 0 to 50 KHz, I could only see low-frequency, close-in spurs, so I zoomed in on a frequency range of 0 to 1 KHz to get a better idea of what these spurs were.

You can see them below.  The spurs are related to 60Hz and its harmonics.  The largest spur, at 120 Hz, is at roughly -90 dBV (dBV means that the measurement is referenced to 1 Vrms).

(Click on image to enlarge)

Below are the actual spur amplitudes from the SigLab plot, above.  The "0 Hz" value represents the DC level driving the EFC input of the 10811 oscillator module.

Will the 60 Hz spurs cause a problem?

Recall from the earlier "Integration Period" table that the 16-bit DAC, to completely cover a span of 10 Volts total, will have a step-size of about 1.53E-4 volts per step.

Our worse spur in the table above is at 120 Hz, and its amplitude is -89.8 dBV (dB relative to 1 Vrms).

These amplitudes are RMS values, so the peak value of the 120 Hz spur is -86.8 dBV, or 4.6E-5 volts.  This is smaller than our minimum DAC step-size of 1.53E-4 volts (by about a factor of 3), so there should be negligible modulation of the oscillator by these spurs.


Phase noise:

Well, EFC looks pretty good, but let's actually measure phase noise.  Here are two measurements taken with my HP 8568B and the KE5FX Phase Noise program.  One measurement is of one of my GPSDO's 10 MHz outputs.  The other measurement is of one of the 5 MHz outputs.

(Click on image to enlarge)

Interestingly, it's showing some spurs around 120 Hz and 240 Hz from the carrier.  If I make a measurement directly on my Spectrum Analyzer, it looks to me that the 120 Hz spur's amplitude is about -80 dBc. 

Here's a screenshot from my 8568B:

(Click on image to enlarge)

The differences in spur-amplitudes between the spectrum analyzer's measurement (above) and the KE5FX phase noise program are probably due to the phase noise measurements being expressed as "dBc per Hz", rather than simply as straight dBc.

Regarding these spurs -- they are much larger while the oscillator's oven is heating up, when significantly more current is being drawn from the power transformer.  So I suspect the transformer's magnetic field is coupling into either the EFC signal, the 10811 internal electronics, or into the output cabling between my MC3487 and the BNC connectors on the back panel (which are routed fairly close to the power transformer).


Other Notes:

Some amount of phase noise will be introduced to the 10811's 10 MHz signal by the logic gates through which the oscillator's signal passes (i.e. the 'LVC1G14, the 'LVC1G80, and the MC3487.  I measured the phase noise of the 10811 oscillator, and it looked to me that its phase noise was a couple of dB less than then phase noise of the GPSDO digital outputs.  For reference, an interesting measurement of  phase noise created by the 74AC family can be found here:  http://www.ke5fx.com/ac.htm


Final Thoughts:

This GPSDO, in my opinion, does a very good job at meeting my needs, but it is a Frequency-Locked loop, not a phase-locked loop, and thus it does not "continuously" adjust the GPSDO's frequency, but instead nudges it one way or another only after a long time interval, sometimes on the order of many hours, during which the frequency will have drifted away from its ideal frequency.

A better approach, if required, would be a phase-locked loop approach (e.g. Brook Shera's design or similar), but this type of architecture will need additional hardware (or so I believe).  The frequency-lock approach (per my design in this blog post) simplifies the hardware design, but at some cost to performance.

Also, I do not know how robust my method of "continuously variable measurement durations" (as described earlier in the "New Algorithm" section) actually is.  You may want to use the VE2ZAZ's fixed-duration approach, instead.


My  GPSDO posts:

(This post) An Arduino-based GPSDO (frequency-locked loop):
http://k6jca.blogspot.com/2016/02/an-arduino-based-gps-disciplined.html

Simulating Brooks Shera Phase-Locked Loop GPSDO (his design that appeared in the July, 1998 issue of QST):
https://k6jca.blogspot.com/2018/12/simulating-brooks-shera-w5ojm-gpsdo.html

Implementing the Brooks Shera Phased-Locked Loop GPSDO on an Arduino Platform:
https://k6jca.blogspot.com/2019/02/an-arduino-version-of-brooks-sheras.html


Thanks!

A big thanks to Dick Benson, W1QG, for allowing me to use his 5370B Matlab plotting tool and for letting me kick around ideas with him.

Standard Caveat:

As always, I might have made a mistake in my equations, assumptions, drawings, or interpretations.  If you see anything you believe to be in error or if anything is confusing, please feel free to contact me or comment below.

Also, if you are modifying the MikroElectronika GNSS2 GPS Receiver board, please remember that you are making these modifications at your own risk. 

And finally, I should add -- this design and any associated information is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.