An FPGA SDR HF Transceiver, Part 1

If you've arrived here, you've arrived at an old post that has been superseded!

To go to the new FPGA SDR Transceiver posts, please click

HERE!

(Thanks!  Jeff, K6JCA)

SDR Notes: The Mixer Mathematics of Digital Down Conversion

(This post originally started as notes to myself, placed here so that I would not lose or forget them.)

The Digital Down Converter (DDC) is a basic block of FPGA-implemented SDR radios.  Per Wikipedia:

A digital down-converter (DDC) converts a digitized, band limited signal to a lower frequency signal at a lower sampling rate in order to simplify the subsequent radio stages. 

Below is an illustration of a DDC stage (don't worry about the math -- it will be explained later in this post):

(Click on image to enlarge)

Note that, for demonstration purposes, the equations assume a simple cosine as the input signal.

This post will look at the math involved in the Mixer.  I will leave the low-pass filter and sample-rate conversion for later posts.


First, some basics...

Basic Down-Conversion Mixing:

Historically, down-conversion has been done in the analog domain with a "real" oscillator signal (e.g. cos(ω_ot)).  Let's first look at this type of basic mixer before tackling mixing in the Digital (and Complex) domain.

For demonstration purposes, I'll simplify the receive-signal to be a just simple cosine signal at a single frequency:  cos(ω_it).

This mixer's math is straightforward.  Let the output of the mixer be y(t).

y(t) = cos(ω_it)*cos(ω_ot)

We can use Trigonometric identities to express this equation in terms of sums and differences of frequencies:

y(t) = (cos((ω_i−ω_o)t))/2 + (cos((ω_i+ω_o)t))/2

But this basic form of mixing has a problem if one would like to down-convert an incoming RF signal directly to base-band:  both the desired signal and its image will appear in the base-band spectrum (which, for the typical shortwave receiver, is the audio spectrum we want to hear).

This type of conversion is known as "Direct Conversion."  For an example, see Rick Campbell's R1 receiver.

Visually, here's how the frequencies shift, shown as a single-sided spectrum (which one would see if one did a "Trigonometric" Fourier Transform on the RF input.  Refer to equation 27 in this paper by Robert Marcus for more detail)::

Note the overlap of the two signals down at base-band...our "desired" signal would be interfered with by any signal at its "image" on the other-side of the oscillator's frequency!

And if we were to consider the Mixer's input to be a band of frequencies:

You can see the spectrum fold-back with the result being the desired signal and its undesired image occupying the same base-band frequencies.  A single-sided spectrum graph has no negative frequencies, and any "calculated" negative frequencies become positive frequencies via the Trig identity: cos(−ω_it) = cos(ω_it).

We can get around this image problem by working in the domain of Complex Numbers (numbers with a real and an imaginary component).  Specifically, we will use exponentials with complex arguments.


Complex Exponentials -- A Quick Review of Some Important Identities:

Euler's Formula:

e^jθ = cos(θ) + jsin(θ)

Similarly,

e^−jθ = cos(θ) - jsin(θ)

From Euler's Formula it can be shown that:

cos(θ) = (e^jθ + e^−jθ)/2

sin(θ) = j(e^−jθ − e^jθ)/2

Let's consider a real signal at a single frequency f_o (Hz), and let's express this signal simply as a cosine: cos(2πf_ot).

Using the identities above this signal can be expressed in exponential form:

cos(2πf_ot) = (e^j2πf_ot + e^−j2πf_ot)/2

Below is a representation, in the complex frequency-domain, of this signal :

(Source:  Quadrature Signals: Complex, but not Complicated, Richard Lyons)

And, for completeness, here are a two more illustrations:

First, to complement the cosine representation above, here is a complex frequency-domain representation of a sine wave signal of the same frequency:

(Source:  Quadrature Signals: Complex, but not Complicated, Richard Lyons)

Second, we can use both of these illustrations of the signals in the complex frequency-domain to generate a visual representation of Euler's Formula, showing how the imaginary sine terms are rotated (when multiplied by j) to the real axis and then added to the real cosine terms:

(Click on image to enlarge)

(Source:  Quadrature Signals: Complex, but not Complicated, Richard Lyons)

Note that multiplying sin(2πf_ot) by j rotates the sine's complex frequency-domain imaginary components "clockwise" by 90 degrees to the real plane. If we were to instead multiply sin(2πf_ot) by -j, the rotation would be 90 degrees counter-clockwise, with the results again being on the real plane, but pointing in the opposite direction, and our final result (after adding the real components together) would represent Euler's Formula for a negative frequency, -f_o, instead of f_o.

Frequency Conversion using Complex Numbers in Exponential Form:

Let's use the above equations to shift signal frequencies.

First, for ease of demonstration, let's define an input RF signal, x(t), to be a simple cosine at a single frequency, f_i Hz.  I.e.:

 x(t) = cos(2πf_it)

To minimize my typing, I will express this signal's frequency as its radian equivalent, where ω_i = 2πf_i.

x(t) = cos(ω_it)

Using the identities above we can express cos(ω_it) as the sum of two complex exponential numbers:

cos(ω_it) = (e^jω_it + e^−jω_it)/2

Now let's now define our oscillator's signal as a complex exponential:

e^jω_ot

(Note that from Euler's Formula e^jω_ot  = cos(ω_ot) + jsin(ω_ot)).

Our next step is to multiply together the input signal and the oscillator signal.  We will call the output y(t).

y(t) =  cos(ω_it) * e^jω_ot

Replacing the input signal with its equivalent complex exponential form, the output equation becomes:

y(t) = [(e^jω_it + e^−jω_it)/2]*e^jω_ot

Finally, let's finish the multiplying.  The result is:

y(t) = (e^{j(ω_i+ω_o)t} + e^{j(−ω_i+ω_o)t})/2

Note that the frequency components in the arguments of both exponentials have been shift up by ω_o.

By the way, we can use Euler's Formula to express this result in terms of sines and cosines:

y(t) = [cos((ω_i+ω_o)t) + jsin((ω_i+ω_o)t)]/2 + [cos((-ω_i+ω_o)t) + jsin((-ω_i+ω_o)t)]/2

You can also derive the same result using only sines, cosines, and trigonometric identities and not complex exponentials, but the math is more cumbersome.


OK, so an NCO's output that is a complex-exponential with an argument whose frequency is +ω_o will shift the incoming signal's frequency up by ω_o.  What happens if the NCO's frequency is negative, i.e. -ω_o?

Using the same approach as above, the DDC's output can be shown to be:

y(t) = (e^{j(ω_i−ω_o)t} + e^{j(−ω_i−ω_o)t})/2

In other words, both of the input signal's complex exponentials are shifted down in frequency by ω_o. Voila, a Digital Mixer!

Again, note that y(t) can be expressed as the sum of sines and cosines.  Using the identities at the start of this post:

y(t) = [cos((ω_i-ω_o)t) + jsin((ω_i-ω_o)t)]/2 + [cos((ω_i+ω_o)t) - jsin((ω_i+ω_o)t)]/2

which is exactly equivalent to y(t) = (e^{j(ω_i−ω_o)t} + e^{j(−ω_i−ω_o)t})/2.


Mixing with Complex Exponentials:

With a view towards understanding SDR building blocks, let's look at a Digital Mixer using the down-conversion math I've just discussed, and which shifts signal frequencies in the digital (rather than analog) realm.

Below is a Digital Mixer block diagram using complex numbers.  Let's drive the RF input with a simple cosine signal, cos(ω_ot):

Note that the oscillator's output is a complex exponential with a negative frequency in its argument.  This negative frequency, as I've discussed earlier in this post, will down-convert incoming RF.

The input signals' spectrum (shown below as two-sided, with negative frequencies, per the frequencies in the complex arguments in the exponential form of the input signal) will shift down:

(Click on image to enlarge)

The illustration below shows a more general case for a band of frequencies at the DDC's input:

(Click on image to enlarge)

Note that, unlike the analog mixer described at the start of this post, this spectrum has both positive and negative frequencies.  In other words, we do not (yet) have an issue with "image wrapping" while we remain in the complex domain.

Another take on the same topic, from "Complex Signal Processing is Not -- Complex," by Ken Martin, showing the same:

(Source: "Complex Signal Processing is Not -- Complex,",Ken Martin)

Note (again) that the input consists of only real, not complex, numbers, but the output Y(t) will be in complex form, with both real and imaginary components.

Regarding the mixer (c) in the image above, we can express its two outputs (y_r(t) and y_q(t)) in terms of sines and cosines.

First, for simplicity of analysis let's define the input x_r(t) to be a simple cosine at frequency ω_r:

x(t) = cos(ω_rt)

Noting that the oscillator frequency is ω_i in this example, the output y_r(t) is therefore:

y_r(t) = cos(ω_rt)*cos(ω_it)

which, using trig identities, expands to:

y_r(t) = cos((ω_r−ω_i)t)/2 + cos((ω_r+ω_i)t)/2

Now let's do the same to derive y_q(t), the imaginary term...noting that −sin(ω_it) = sin(−ω_it), y_q(t) becomes:

y_q(t) = sin((ω_r−ω_i)t)/2 − sin((ω_r+ω_i)t)/2

Is this result equivalent to the result derived (above) using complex-exponential equations?

Essentially, yes it is!

Let's multiply the y_q(t) channel by j (i.e. add a phase shift of 90 degrees with, say, a Hilbert Transform) and call this output y_q(t)' (note the 'prime').  Its equation is:

y_q(t)' = jsin((ω_r−ω_i)t)/2 − jsin((ω_r+ω_i)t)/2

If we then add y_r(t) and y_q(t)', the result is:

y_r(t) + y_q(t)' = [cos((ω_r−ω_i)t) + jsin((ω_r−ω_i)t)]/2 + [cos((ω_r+ω_i)t) − jsin((ω_r+ω_i)t)]/2

Which, from Euler's Formula, we know is exactly equivalent to:

y_r(t) + y_q(t)' = (e^{j(ω_r−ω_i)t} + e^{j(−ω_r−ω_i)t})/2.

and which is exactly equivalent in form to the equation we derived earlier:

y(t) = (e^{j(ω_i−ω_o)t} + e^{j(−ω_i−ω_o)t})/2

In other words, the trigonometric signals representing y_r(t) and y_q(t) are the fundamental mathematical elements (when the sine terms are multiplied by j) that constitute the Digital Mixer's operation when expressed in complex-exponential form.

Links and Other Resources:

An Excellent Tutorial on Complex Numbers:  Quadrature Signals: Complex, but not Complicated, Richard Lyons

An interesting site to surf:
https://www.dsprelated.com/freebooks/mdft/Positive_Negative_Frequencies.html

Representing a signal as the sum of two exponentials:  http://ece.mst.edu/media/academic/ece/documents/classexp/ee216labs/EE216_Lab4.pdf

The Significance of Negative Frequencies in Spectrum Analysis, Robert Marcus, IEEE Transactions on Electromagnetic Compatibility, Dec. 1967.  (A copy can be found here).

Good tutorial on Complex Number in Radio applications:  Complex Signal Processing is Not -- Complex, Ken Martin

Standard Caveat:

I could have easily made a mistake in any of the above text, equations, or drawings.  If something looks wrong or is confusing, please feel free to contact me.  Thanks!

More Notes on Directional Couplers for HF -- The Twin-Lead "Twin-Lamp" SWR Indicator

The Twin-Lead "Twin-Lamp," by W4HVV

At a recent swapmeet someone asked me how the venerable dual light-bulb "twin-lead" SWR indicator worked.  Although I was familiar with the device, I hadn't looked into its theory of operation, and so I said that it was probably similar to the Monimatch (described here) because, even though it used a length of twin-lead as its sensing element, that length was very short compared to an HF wavelength and thus it could be analyzed from a "lumped-element" perspective using capacitive and inductive coupling, rather than as a transmission line.

But the question made me curious as to what the analysis would be, and I decided to delve deeper...

The 'Twin-Lead "Twin-Lamp"' SWR indicator was a very simple device designed for use with twin-lead (or open-wire) transmission lines, and it used two inexpensive flashlight bulbs as Forward and Reflected Power indicators.  To my knowledge, it was first mentioned in an article by Charles Wright, W4HVV, in the October, 1947 issue of QST (The "Twin-Lamp" -- you can find a copy of the article here).  It then appeared in later editions of the ARRL's "The Radio Amateur's Handbook (up through 1961), as well as in editions of the "Radio Handbook" published by Editors and Engineers (such as the 13th edition, published in 1951, which I have).

The W4HVV's "Twin-Lamp" SWR Indicator is pictured above.  Below is a sketch of a similar design, for the 13th Edition of the Radio Handbook, published by Editors and Engineers:

In other words, the SWR indicator consists of a short length of twin-lead with two lamps attached to it.  The lamps can be attached at either end, or in the center.

Per W4HVV, the configuration shown above as (A) has greater sensitivity than configuration (B), no doubt due to the direct (rather than capacitive) connection of the top of the detection-loop to the transmission line.

W4HVV does an excellent job explaining the detector's operation, and I've included it below (from the original QST article, via the RF Cafe website).  First, the figures he references:

(Click on image to enlarge)

Figure 1

And here is his explanation (I've added some formatting changes for clarity):

Referring to Fig. 1-A, a current, I_L, in the line would induce a current, I₁, in a loop near the line, as shown. If the reactance of the loop is small compared to the resistance of the bulbs A and B, the current I₁ will lag I_L by 90°. This current will, of course, be the same through lamps A and B, and will cause them to burn with equal brightness if they are identical.

Now from Fig. 1-B, we see that bulbs A and B are across the line and in series with a small capacity C. This capacity is, of course, the distributed capacity between the loop and the line. If the reactance of this capacity is large compared to that of A and B the current I₂ will flow and will lead the voltage across the line by 90°. If A and B are identical the current will divide equally between them.

Since I₁ lags I_L. by 90° and I₂ leads E_L by 90°, it is apparent that if I_L and E_L are in phase with each other, I₁ and I₂ will be exactly out of phase.

Fig. 1-C is a combination of the circuits explained above. Condenser C is the capacity between the wires of the loop and the line. Currents I₁ and I₂ are shown as they appear in Figs. 1-A and 1-B. It is now evident that bulb A will light from the sum of I₁ and I₂ and bulb B will light from the difference between these two currents. This is the case for a wave traveling toward the right. In the case of a wave traveling toward the left, the currents will add in bulb B and tend to cancel in bulb A. Thus the device is a form of "directional coupler." When the line is terminated on the right-hand side (marked "load" in Fig. 1-C) by a resistance equal to the characteristic impedance of the line, there is no reflected wave and only bulb A will light. If the load is something different, there will be some reflected energy, and lamp B will burn along with A, the relative brilliance depending upon the relative magnitudes of the transmitted and reflected energy. These facts are what make the device so useful as a standing-wave indicator.

In the foregoing discussion, three conditions were set up: 

1. bulbs A and B should be identical; 
2. the reactance of the loop should be small compared to the impedance of A and B; and 
3. the reactance of the coupling capacity should be high compared to the bulb impedance. 

To satisfy the first, bulbs of the same characteristics were used, and in the interests of sensitivity, these were 2-volt 60-ma. flashlight bulbs. For the second and third considerations, the length of the coupling loop must be kept short compared to a wavelength. It was found that, for 50-Mc. operation and a transferred power of about 20 watts, a loop length of about 4 inches was a good compromise between sensitivity and the satisfaction of the above conditions. For 28 Mc. it can run a few inches longer, and at 144 Mc. an inch or so shorter. In any event, the length is not critical.

A Closer Look...

Let's dive a bit deeper and see if we can't derive some equations to describe the operation of the "Twin-Lamp".

I'm going to use the principle of Superposition (a fundamental principle in Linear Network Theory and an essential tool for circuit analysis).  I'll first look at how the voltage across the transmission line creates currents in the two lamps of the detector (independent of the current through the transmission line).  I'll then look at how the current through the transmission line (independent of the voltage across the line) creates its own currents in the two lamps.  And then I'll combine these two currents, calculated independently, to get the final aggregate lamp currents.

First, let's consider a twin-lead (or open-wire) transmission line driven by a source at one end and terminated by a load (Z_load) at the other end:

If Z_load differs from the characteristic impedance of the transmission line (Z_o), then there are reflections along the transmission line and, if we were to measure current and voltage at an arbitrary point along the transmission line, we would see an impedance, Z_load' (i.e. Z_load prime), where the value of Z_load' changes as we move along the transmission line, and which only equals Z_load at distances equal to integer half-wavelengths as we move away from the load back towards the generator.

(Click on images to enlarge)

So, if we were to measure V and I at the (arbitrary) point shown above, it would appear to us that we had a load of impedance Z_load' attached at that point, so that the transmission-line circuit can be modeled now as a much shorter line with a new load, Z_load':

where:  Z_load' = V/I, measured at that point.

What happens if we put our "Twin-Lamp" SWR indicator at that point, so that it measures the same V and I that creates Z_load' ?  (Remember, the length of the Twin-Lamp detector should be significantly shorter than the wavelength of the signal being measured to ensure that the I and V it measures are essentially the same as those at Z_load' ).

(Click on image to enlarge)

First, let's look at the currents through the lamps due to the voltage across the transmission line (I_lamp-C).  We can consider the bulbs to be capacitively coupled to the transmission line and, if the bulbs have identical impedances, then the currents through them will be equal and in the same direction:

(Click on image to enlarge)

The circuit is a simple voltage divider, and thus the current I_c is equal to the voltage across the line at that point divided by the sum of the series-impedances across that point:

I_c = V/(1/jωC + (Rlamp || Rlamp)+ 1/jωC)

I_c = V/(Rlamp/2 + 2/jωC)

if (2/jωC) >> (Rlamp/2) -- per W4HVV's third condition, above (that is, if the capacitive reactance is much greater than the lamp resistance), then the equation reduces to:

I_c = V*jωC/2

Assuming the impedances of the two lamps are identical, then I_c splits equally through each lamp branch:

I_lamp-C = I_c/2 = V*jωC/4

Now let's look at our second lamp-current factor: lamp current due to transmission-line current (I_lamp-L).  The transmission-line current inductively couples to the detector loop:

(Click on image to enlarge)

If the reactance of the mutual coupling is very small compared to the bulb resistances (per W4HVV's second condition, above), then this coupling can be modeled as two voltage sources (refer to this post as to why):

(Click on image to enlarge)

The lamp current due to inductive coupling, I_lamp-L is simply the current in the detector's loop, and it can be derived as follows:

I_lamp-L = V_induced/(R_lamp +  R_lamp)

therefore,

I_lamp-L = 2*jωM*I/(2*Rlamp)

substituting V/Z_load' for I and reducing,

I_lamp-L = jωM*V/(Rlamp*Z_load')

So, we've calculated the currents due to capacitive coupling and due to inductive coupling.  Now let's combine them and calculate the voltage across each lamp:

(Click on image to enlarge)

Notice that the two current-components through lamp-F are in the same direction.  Therefore they add, and the voltage across this lamp can be calculated as:

V_lamp-F = (I_lamp-C + I_lamp-L)*R_lamp  

But the current-components through lamp-R are in opposite directions and therefore subtract.  The voltage across this lamp can be calculated as:

V_lamp-R = (I_lamp-C - I_lamp-L)*R_lamp     

If the transmission-line is properly terminated in a load equal to the line's characteristic impedance Z_o, then the SWR should be 1:1 and lamp-R should be extinguished -- the voltage across it should be zero.  Therefore, substituting Z_o for Z_load' and setting V_lamp-R to 0 volts, our last equation becomes:

0 = (I_lamp-C - I_lamp-L)*R_lamp 

    

0 = (V*jωC/4 - jωM*V/(Rlamp*Z_o))*R_lamp 

Which reduces to:

0 = C/4 - M/(R_lamp*Z_o) 

This last equation gives us the relationship we need to design the appropriate values of capacitive and inductive coupling between the "Twin-Lamp" detector and the transmission line:

4*M/C = R_lamp*Z_o 

Note the similarity between this equation and the equation derived for the monimatch:

M/C = R1*Z_o 

Other notes:

The analysis above (and W4HVV's conditions) assume that the resistance of the two bulbs is identical.  But a bulb's resistance should increase as its filament heats up, so there will be a difference in the resistances of an illuminated "Forward" lamp compared to a dark "Reflected" lamp.  This difference shouldn't affect the capacitive coupling currents (as long as the capacitive reactance is significantly greater than the "hot-filament" resistance), but it will affect the current due to inductive coupling.  In this case (of unequal lamp resistances), I_lamp-L becomes:

I_lamp-L = 2*jωM*V/((R_lamp-F + R_lamp-R)*Z_load')

and thus the final equation changes to:

8*M/C = (R_lamp-F + R_lamp-R )*Z_o 

How will this difference affect the actual SWR reading?  I don't know.

As for the other conditions mentioned by W4HVV, if this device's length is kept short, then the overall capacitive coupling should be fairly small, and W4HVV's third condition of capacitive reactance being significantly larger than lamp resistance should be met.

As for W4HVV's second condition, that loop-reactance should be kept small compared to the reactance of the bulbs, I'm not sure how well this condition is met in this design, but as long as the "Reflected" bulb is dark when the line is terminated with a resistance equal to the line's characteristic impedance, Z_o, I wouldn't worry too much -- after all, this device is not meant to be a laboratory-grade instrument!

And when verifying the operation of this type of SWR Indicator by terminating the transmission line with a resistive (non-inductive) load whose value equals the transmission line's characteristic impedance, Z_o, (to ensure that the "Reflected" bulb is dark under "match" conditions), note that the actual Z_o of twin-lead or ladder-line might differ from its assumed value.  For example, 450-ohm ladder-line might be closer to 400 ohms.

But, as W4HVV points out, his "Reflected" bulb typically did not begin to light up until SWR was at least 1.5:1, which represents a significantly larger difference in impedances than 400 ohms compared to 450 ohms.

OK, that's it for this analysis!


Some Twin-Lead SWR Indicator articles:

The "Twin-Lamp," Wright, W4HVV, QST, Oct., 1947

"The Eyes Have It," Paddon, VE3QV, QST, Oct., 1948

"The Coax Twin-Lamp," Keay, W0SJK, QST, Nov., 1948

"An Improved Twin-Lamp", Fisher, VE3ALQ, QST, Oct., 1949

"Measuring Center Impedance of Antennas with the "Twin-Lamp," Gross, W2OXR, QST, May, 1950

"The Balanced Twin-Lamp," Wood, K2BUZ, QST, Nov., 1956

"A Reflectometer for Twin-Lead," Brown, W6HPH, QST, Oct., 1980



Links to my Directional Coupler blog posts:

Notes on the Bruene Coupler, Part 2

Notes on the Bruene Coupler, Part 1

Notes on HF Directional Couplers (Tandem Match)

Building an HF Directional Coupler

Notes on the Bird Wattmeter

Notes on the Monimatch

Notes on the Twin-lead "Twin-Lamp" SWR Indicator

Calculating Flux Density in Tandem-Match Transformers


And some related links from my Auto-Tuner and my HF PA posts:

Auto Tuner, Part 5:  Directional Coupler Design

Auto Tuner, Part 6:  Notes on Match Detection

Auto Tuner, Part 8:  The Build, Phase 2 (Integration of Match Detection)

HF PA, Part 5: T/R Switching and Output Directional Coupler


And other links:

Modeling RF transformers


Final Caveats:

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

Use a Dremel Router to Prototype Circuits with SMD Parts

(Little pads for a little part)

My projects are almost always "one-off" projects.  That is, I build one and that's it.

And, although I've overseen the layout of many PCBs for my designs over my career as an engineer, I never make a PCB for my own personal projects.  I find that it is much faster for me to just build the circuitry on a piece of "vector-board" (i.e. a board with or without plated-through holes, usually on 0.1" centers) or "dead-bug" style on a piece of unetched copper-clad PCB material.

(Another construction technique is known as "Manhattan-style", which requires gluing (or soldering) small pieces of PCB to an underlying copper-clad PCB base.  I don't use this technique because, to me, it has always seemed to be a bit too much work, and I have also heard that, if gluing the pcb "pads" to the base, heat from a soldering iron can sometimes cause the glue to soften and the pad to shift.)

For me, my "vector-board" and "dead-bug" techniques work great for most of my projects, unless those projects use surface-mount devices.

SMD resistors and caps aren't usually a problem, as they can be easily soldered between two pads on a vector board.  But there can be issues when trying to incorporate multi-pinned SMD ICs.  I find that the tiny pin-to-pin spacing is difficult to work with, especially when trying to solder wires or components to them.

But SMD has much to offer the circuit builder, not least of which includes smaller circuits and better EMI and Susceptibility performance.

I've tried creating pads in the PCB copper using an Xacto knife to cut through the copper and then removing the unwanted metal with the hot tip of a soldering iron.  But this is time consuming -- I need to make several passes with the knife before I've cut deep enough into the copper to allow its removal with the soldering iron's tip.

But there's a better way...

Recently I was thumbing through a catalog of luthier tools from Stewart MacDonald and I came across a Precision Router Base for a Dremel tool:

With the additional mass and stability it would add to a Dremel tool, this base looked to me to be perfect for cutting SMD pads on PCB material.

(Click on image to enlarge)

So I ordered one, gave it a try, and was happily surprised at the results.

Here's what I do...

The Procedure:

1.  I first draw the outlines of the pads (or traces) on the PCB stock:

2.  Then, before routing, I verify that the pads (or traces) have the correct locations and sizes:

3.  After mounting the bit in the Dremel tool and adjusting the router-base height so that the bit just barely protrudes below the base, I then carefully route along the outlines I've drawn (in the example shown below I'm using a 0.1mm, 30 degree Engraving Bit with the Dremel tool).  Note there is some wiggle on some of the lines -- they aren't perfect.  But these cuts are much straighter than if I had tried to "freehand" the routing without using this base.  (And they could probably be improved upon if I used a straight-edge to guide my router-base):

4.  I can then mount the part (i.e. solder it to the pads).  However, before soldering, first verify with an ohmmeter that all pads are isolated from one another:

That's all there is to it.  Pretty simple! 

What bit to use?

I experimented with four different bits:

1. 1/32" Diameter Straight Router Bit (Stewart MacDonald item 1187)
2. Tapered Router Bit (Stewart MacDonald item 1180)
3. 30 degree, 0.1mm Tip Diameter Conical Engravers Cutter V-bit.  (Note: you can find these on Amazon and eBay for less than the Endpoint version identified in the link above).
4. Dremel 105 Engraving Cutter Bit (1/32").

Of these four bits, the "30 degree 0.1mm" tip gave me the finest engraved lines, making it best suited for cutting SMD pads.

The Dremel 105 "Engraving Cutter" bit (which has a cutting-ball at its tip end) also works well, but it results in a broader line -- useful if you need more pad-to-pad isolation (i.e. if creepage and clearance distances are issues).  Note, too, that there are other "Engraving Cutter" bits with larger diameter balls at their ends, should you desire even greater isolation between pads or traces.

The two Stewart MacDonald bits were less than satisfactory.  The first reason is that they require a 3/32" Dremel collet (which I discovered I did not have, and so I needed to purchase one).  These two bits result in rougher edges along the cut (in fact, I took two passes for the 1180 bit route in the photo below -- it was much rougher with only one pass).  Plus, being flat-bottomed, a cut is a little bit harder to start with the Stew-Mac router bits, compared to the 105 (ball-tip) bit or the 30 degree, 0.1mm tip.

(Click on image to enlarge)

Other accessories:

I also purchased a companion Mini Air Pump from Stewart MacDonald to blow away the dust as I was cutting the pads.  It helps, somewhat, but a significant amount of debris is still left on the board and it can obscure what I'm trying to cut, so the pump's usefulness is limited.  A more powerful air pump would probably help.


That concludes this post.  If this construction technique interests you, I would recommend experimenting with different tips and Dremel motor-speeds to see which combination works best for you.


Other Notes:

This method works well with SOIC pin spacing.   But I personally find TSSOP pin spacing to be a bit too tight to use this method.


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, when using a Dremel tool, be sure to always wear protective eye-wear.  This caution is especially important when using some of the fine-tips I mention above, as it is possible to break them if one isn't careful, and you do not want a piece of a tip flung into your eye at high speed.

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.

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.