Simulating the Brooks Shera (W5OJM) GPSDO Algorithm

Note, 6 February 2019:  I've implemented the Brooks Shera PLL algorithm in my Arduino GPSDO.  See this blog post:  http://k6jca.blogspot.com/2019/02/an-arduino-version-of-brooks-sheras.html

Introduction:

The July, 1998 issue of QST had an interesting article by Brooks Shera, W5OJM, titled, "A GPS-Based Frequency Standard."

I thought this GPSDO (GPS Disciplined Oscillator) would be great addition to my home lab, and so I built his design and used it to lock my venerable HP-106B Frequency Standard to GPS:

Many years later, I decided to build a second GPSDO for my "home away from home" lab.  This one, though, was a simpler design, and it used a frequency-locked loop in lieu of a phase-locked loop to discipline the internal oscillator:

(I describe this GPSDO in the blog post:  http://k6jca.blogspot.com/2016/02/an-arduino-based-gps-disciplined.html)

This frequency-locked GPSDO works quite well for my application (its accuracy is typically better than 0.1 parts-per-billion  -- i.e. better than 0.1 Hz at 1 GHz).  But I got to wondering -- what would it take to modify my Frequency-Locked Loop GPSDO into a Phase-Locked Loop GPSDO?

My knowledge of control-loops in general and phase-locked loops in particular was sketchy.  So it was time to do some research!

An obvious starting point was to investigate Brooks Shera's original design.  I had the PIC source code for a version of his GPSDO software, but this source code is only in the form of a print-out -- the file from which it had been printed had been lost long ago.  And, unfortunately, Brooks is now a Silent Key.

The source code holds the keys to understanding the algorithm.  Unfortunately, it's in PIC assembly language!

(Shera's code for his IIR filter)

I know nothing about PIC assembly language coding.  But with Google at my fingertips and Brooks' copious comments in the code itself, it wasn't too difficult to work out what was going on.


Simulink PLL Phase-Domain Modeling:

The Mathworks Simulink program is a great tool for modeling and understanding systems, and it is this program that I used to simulate Brooks Shera's design.

With advice and help from W1QG, I modeled Shera's PLL in the "Phase Domain".  An example of this type of modeling is shown, below:

You can find more on Phase-Domain Modeling on the Mathworks' website:
https://www.mathworks.com/company/newsletters/articles/modeling-and-simulating-an-all-digital-phase-locked-loop.html

Simulink Model of Shera's design:

I created the Simulink model, below, from Shera's PIC code:

If you would like to play with this model, you can download it from GitHub:
https://github.com/k6jca/Brooks-Shera-GPSDO.

By the way, the different colors refer to different timing domains.  This model has its timing expressed in terms of "samples" (rather than seconds).  The Red paths and blocks operate every sample (which, in real life, is at a 1 sample-per-second rate), while the Green blocks operate at a sample-rate that is 1/30th the "Red" rate (i.e. one sample every thirty seconds).

The color-code is described in this legend:

Let's take a closer look at the functional subsections of this model and Simulink blocks that comprise them.

Modeling the Phase Detector:

Overview:

The Phase Detector subsection consists of the following parts of the Simulink model:

In Brooks Shera's GPS Disciplined Oscillator, the Phase Reference is the GPS Receiver's 1 PPS (1 Hz) pulse.  This reference is represented by the "Phase Reference" block in the upper left-hand corner of the simulation model, above.

In real life, the GPS 1 PPS pulse is jittery.  So, to reduce the effect of this jitter, Shera's design, over a period of 30 seconds, sums together, second-by-second, each second's phase-delta.  After 30 seconds of summing, this "aggregated" sample (representing  the sum of 30 seconds of data) is then the sample that is sent to his Low-Pass Filter block.

Therefore, once every 30 seconds, an "aggregated" sample is applied to the input of the low-pass filter, whose output is then sent, once every thirty seconds to the DAC which drives the VCO.  In other words, the system operates at a one-sample-every-thirty-seconds sample rate.

To measure phase delta for a 1 second period, Shera's circuit, in hardware, counts the number of 24 MHz clocks that occur during the period of a gating pulse that is set at the rising edge of the GPS receiver's 1 PPS pulse and reset at the rising edge of a 312.5 KHz signal derived from the GPSDO's VCO output (i.e. 10 MHz / 32).

Here is my representation of his Phase Detector circuit:

The duration of this gating-pulses' period can range from 0 to 3.2 microseconds, depending upon when the rising edge of the 312.5 KHz signal occurs in relation to the 1 PPS rising edge.

At its longest duration of 3.2 us, the gating pulse would allow almost 77 (76.8) of the 24 MHz clocks to be counted.  Therefore, over 30 seconds, if the duration of this gating pulse were to continue at its maximum of 3.2 us,  a total of 2304 clocks would be counted.  Thus, 2304 represents the MAXIMUM phase delta value from the phase detector.  And zero is its minimum phase delta value (when the rising edges of both the 1 PPS and the 312.5 KHz signals are co-incident).

In the Simulink model the phase detector's output (from the Phase-wrapper) ranges from 0 to 1.  So, to simulate Shera's second-by-second phase-delta count, I multiply the Phase-wrapper's output by 76.8.

These second-by-second phase delta samples are then "aggregated" by a 30-sample FIR filter whose coefficients are all set to 1.    Therefore, if there are 30 consecutive single-second phase-delta samples that are each equal to 76.8, the output of the FIR filter will by 2304.

This phase-delta is compared against a "setpoint".  The GPSDO's control loop will attempt to adjust the VCO frequency so that the phase output from the phase detector equals the setpoint -- that is, so that this difference between the two is zero.  I've chosen the setpoint to be mid-range between the phase-detector's minimum and maximum values (i.e. set to 1152) -- this will give the most "headroom" for the phase to move around (without wrapping around) if there is a severe perturbation to the system.

Note that Brook Shera's code uses a setpoint of 800 for revision 1.29 and later, and 1024 for earlier revisions.  I simply set mine to 50 percent of the range (1152).

Let's look more closely at some of the blocks in the diagram, above.

Perturbation Block:

The Perturbation block allows me to add a perturbation the reference phase at a specific time, so that I can see how the system responds to it.

Clicking on the "Perturbation" Block at the model's top level brings up the following window:

The range of values are from 0 to 1 (representing 0 to 100% of the Phase Delta's range).  Note that the model's phase reference is set to 50% of its range (i.e. actual value is 0.5), and this Perturbation block adds to it (at sample-step 100) a step-function of amplitude 0.49.

Phase-wrapper Block:

The Phase-wrapper block essentially "wraps" the phase around if the phase were to be greater than 360 degrees (i.e. 2 pi radians) or less than 0 degrees.  It is instantiated as a MATLAB function, using Simulink's "MATLAB Function" block.

Clicking on this block pulls up its (function) code:

The output from the Phase-wrapper is normalized to 1.  It will then be scaled so that the phase delta equals 76.8 when the Phase-wrapper's output is at its maximum value of 1.

FIR Filter:

Shera's design then aggregates 30 seconds of phase-delta samples.  The Simulink model does this using a simple 30-sample "brickwall" FIR filter whose coefficients are all equal to 1.

Here's the definition of the FIR block:

Down-sampler:

Following sample-summing, Shera's PLL operates at a sample rate of one sample every 30 seconds.

The Simulink model accomplishes this rate conversion with a Down-sampler block, which samples the output of the FIR filter once every 30 samples..

Clicking on the Down-sampler block pulls up its parameter window:

"D" is the amount by which the incoming data stream is down-sampled by.  In this case, it is set to 30, but to use it globally in the model, I've called it D and defined D to be 30 in the model's "InitFcn" callback (under the "Callbacks" tab in the " Model Properties" window):

Modeling the Low-Pass Filter:

Shera allows the user to select one of seven different low-pass filters to filter the phase-delta error signal.

The first filter is what he calls a "Type 1" filter, and it is simply the phase-delta error multiplied by a constant (in this case, 32).

The six other filters are all IIR (Infinite-Impulse Response) low-pass filters whose time constants successively increase by a factor of 2.

Let's look at the IIR filter block...

IIR Filter:

The IIR Filter is the heart of Brooks Shera's GPSDO algorithm.  Here is the Simulink version of his IIR filter as I've gleaned it from his PIC source code:

This IIR filter implements the following equation:

o(n) = o(n-1) + i(n)*(1/F1 + 1/F2) + i(n-1)*(1/F1 - 1/F2)

Where:

• i(n) = input sample (phase-delta) at sample-time 'n'
• i(n-1) = previous input sample
• o(n) = output sample at sample-time 'n'
• o(n-1) = previous output sample
• F1 is a constant, defined by Shera to be between 2^11 and 2^16.
• F2 is a constant, defined by Shera to be 64.

The filter's output, o(n), is then multiplied by the constant 'Kcpu'.  (Kcpu is defined by Shera to be between 2^10 and 2^5, depending upon filter time-constant).

Thus, the DAC input is:

DAC_IN(n) = Kcpu*[o(n-1) + i(n)*(1/F1 + 1/F2) + i(n-1)*(1/F1 - 1/F2)]

Modeling the DAC and the Controlled Oscillator (VCO)

The output of Shera's low-pass filter stage is the input to his 18-bit DAC, which generates the voltage to control the VCO's frequency, as shown below.

Shera has defined for his algorithm a VCO "sensitivity" S (delta-F/F/Volt, as referenced to the voltage at the DAC output pin) that must be 7.5E-9 per volt (i.e. 7.5 ppb per volt).

This means that a 10 MHz oscillator should have a Kv equal to 0.075 Hz/V.  (Kv = S*VCO_Frequency = 7.5E-9*10E6 = 0.075 Hz/V).

Given the wide range of Kv's that VCOs could have, Shera's design has a resistive attenuator between the DAC's output and the VCO's EFC (Electronic Frequency Control) pin that is used to reduce a VCO's higher Kv to provide the target sensitivity of 7.5 ppb per volt.

For example, my HP VCO module (p/n 10811-60111) has a measured Kv of 0.32 Hz/volt.  To bring this down to the target Kv of 0.075 Hz per volt, the EFC voltage from the DAC should be reduced by a factor of 4.27.  (Or, you can think of this as multiplying the DAC output by a factor equal to 0.075/0.32, or 1/4.27).

DAC Model:

I've modeled Shera's 18-bit DAC with the following block:

This model simply converts the input (limited to the range of -131,072 to +131,071) to a value between -3 and +3 (representing the output of the Burr-Brown PCM-61P 18-bit DAC).

Note that the limiter is implemented with a MATLAB function block.  Also, I don't bother to quantize the input to be an 18-bit integer.

VCO Model:

The heart of the VCO is an integrator which converts VCO frequency into phase for this "Phase Domain" GPSDO model (note that phase is the integral of frequency):

Finally (and not shown in the VCO model, above) , the VCO output is divided by N to create a 312.5 KHz signal that feeds back to the Phase Detector.  In this Simulink model, the VCO frequency is 10 MHz, and so N is set to 32.


Simulation Results:

Let's see how the PLL reacts to a step-function perturbation of 0.49 on the "phase reference" (i.e. 1 PPS) signal:

Let's first select the "Type 1" Filter (by clicking on the model's switch at the output of the IIR filter block).

Here is the VCO output in response to an input phase perturbation that is a step-function of amplitude 0.49 (i.e. 49% of the maximum phase delta possible):

And here is the DAC input, for the same step-function input:

Note that the peak DAC input is just a bit more than one quarter of full-scale (the latter being 131,071).

There are six possible IIR filters.  Here is the response of the filter to a step-function of 0.49 amplitude with the fastest transient response (i.e. F1 = 2^11, F2 = 64, Kcpu = 2^10).  (All other IIR filters have slower transient responses).

Note that there is some overshoot and a slight amount of ringing.

Here's the DAC input given the same step function and IIR filter:

If we now select the next IIR filter (F1 = 2^12, F2 = 64, Kcpu = 2^9), the response to the step-function input is the same, but stretched out in time by a factor of two:

The DAC input for the same stimulus is:

You can see that the peak DAC input for this latter example is half that of the IIR whose parameters are:  F1 = 2^11, F2 = 64, Kcpu = 2^10.

With each successive filter, Shera increases F1 by a factor of 2 and decreases Kcpu by a factor of two.  The result is that each successive filter has its step-response stretched in time by a factor of 2 from the previous step-response, and its maximum DAC input is halved.


Other Notes and Comments:

1.  Simulation Experiments:

The model allows the user to experiment with the design by changing model parameters.

For example, as a baseline let's again look at the step-response of Shera's first "stock" Type 2 IIR filter (F1 = 2048, F2 = 64, Kcpu = 1024):

If I increase Kcpu by a factor of 2 but keep the other parameters unchanged, the response has less ringing and overshoot, but the rise-time is faster:

I can lengthen the rise and settling times by a factor of two by increasing F1 by a factor of 2 and decreasing Kcpu by a factor of 2 (just like I did earlier, for the "stock" Shera filters), as you can see in the plot, below (I've increase F1 from 2048 to 4096 and I've decreased Kcpu from 2048 to 1024).

Note that the response still has the same shape as the previous response.  It's just "stretched out" by a factor of 2, indicating that it has a lower cutoff frequency.

2.  Using Simulink to perform Control System Analysis:

For the plots, above, I used a perturbation of 0.49, manually created by me, that was added to the reference phase.  But it is possible to have Simulink make various types of linear system analyses (e.g. step-response, impulse response, Bode plot, etc.) without the user actually having to create a stimulus such as I did.

To do this, you first need to insert into the model the test points where Simulink will, for example, insert an input perturbation and where it will measure the response.

Note: to use this capability of Simulink, I believe you must purchase the MathWorks Simulink Control Design product (other toolboxes of theirs might be required, too).

To insert these test points, the signal line into which each test point is to be added must first be high-lighted and then right-clicked upon to select the menu to then further select "Linear Analysis Points".

In my model, I'm using the two points for signal-insertion and measurement shown below:

I can then select which type of Linear Analysis I'd like to perform by going through the following menus...

...and then selecting "Linear Analysis..." at the last menu, which opens the window for the Linear Analysis Tool:

You can see along the top of this menu some of the types of analysis that one can do.  For example, if I click on the Bode plot, the Linear Analysis Tool generates the following Bode plot.

Clicking on Step Response produces the following plot:

In shape, it looks like the one I generated with my manually-inserted step perturbation.

I'll change the IIR filter to the filter that is the next step up (i.e. lowering the cutoff frequency by a factor of 2) and generate a new Bode plot:

And yes, the cutoff frequency has shifted down by a factor of two from the previous Bode plot.

Here's the Bode plot for the Type 1 filter:

And the Type 1 filter's step-response:

3.  Adding a jittery 1 PPS reference to the model (Quectel L76 GPS Receiver):

My second GPSDO with the Frequency-locked loop (FLL) uses a Quectel L76 GPS Receiver.  This receiver's long-term variation in jitter seems to be worse to me than the jitter on a Trimble Resolution-T receiver.

I was curious how Shera's design would handle the actual jitter on the Quectel receiver.  But first, to simulate its effects, I needed capture actual jitter and put it in a file for Simulink to use.

I built a simple phase-detector for my FLL GPSDO whose output could be read by the Arduino processor.  Essentially, the phase-detector charges a capacitor between the rising edge of the 1 PPS signal and the rising edge of 2.5 MHz derived from my VCO (10 MHz / 4).  So the interval between these two edges can range from 0 ns to 400 ns.

Here is my Phase Detector circuit:

Similar RC charging techniques are used here (Lars GPSDO) and here (Pedersen).  (Note that if I were to incorporate this circuit permanently into my GPSDO, I would probably add one more divide-by-two stage to create 1.25 MHz and thus lengthen the maximum Phase Detector period from 400 ns to 800 ns, to get more "headroom" for 1 PPS edge variations.)

After the capacitor has charged, its voltage is read by the Arduino's ADC.  For my chosen RC time constant and ADC reference voltage (1.1V), the maximum count from the ADC is about 800 for the maximum period of 400 ns.  In other words, each count increment represents a time-interval of 0.5 ns (400ns/800counts).

To capture jitter data, the FLL GPSDO's DAC output is fixed (held), so that the only variation in phase jitter (ideally) comes from the GPS 1 PPS pulse.

The Arduino's ADC data is sent, second by second, to a serial port, where I capture it and place it into an Excel spreadsheet.  I can then manipulate the data -- for example, converting ADC counts to nanoseconds of jitter.  Once transformed into a nanosecond format, it can be used as a perturbation source for my Simulink simulation.

Note that, although the DAC output was fixed, there is some drift because its output voltage was not set to produce exactly 10 MHz from the VCO.

To bring this Excel data into the model, I placed a "From Spreadsheet" block and identified in its parameters which spreadsheet (and data therein) to use.

Here is the portion of the Simulink model that I modified.  I've circled in red the "From Spreadsheet" Simulink block:

Clicking on the "From Spreadsheet" block opens up its block-parameters window:

You can see that the data comes from Sheet 3 of the specified .XLSX file, and that the first column of the data is "Time".

Here's the jitter data (in nanoseconds) in the spreadsheet:

(I manually zeroed-out the first 30 values of the normalized phase-delta to start the sequence with no noise).

This jitter data, in nanoseconds, is then normalized in the Simulink model to Shera's phase-detector period of 3200 ns.  I added a "manual" switch to the Simulink model to let me choose the Excel noise file or my original step-function perturbation.

Here's how the .XLSX data looks after it has been normalized and summed with the 0.5 "reference" phase (but before it has been multiplied by 2pi to drive the Phase Detector -- i.e. looking at the input to the 2*pi gain block).  Note that, being normalized, the limits of its range are 0 to 1, and that it should normally sit around 0.5, chosen because this is the VCO's target phase, given the setpoint value of 1152.

Let's take a look at Shera's VCO response to this data...

Here's the VCO Response with Shera's Type-1 Filter selected:

And here's the VCO response with the first IIR filter selected (F1 = 2048, Kcpu = 1024):

Some cosmetic differences, but not overwhelming.

Let's try a more severe IIR filter.  Here's the fourth filter of the six (F1 = 16384, Kcpu = 128).

Quite a bit different!

4.  Adding a jittery 1 PPS reference to the model (Trimble Resolution T GPS Receiver):

I recently purchased a Trimble Resolution T GPS receiver from eBay, and I thought I would compare its 1 PPS jitter to my Quectel L76's 1 PPS jitter.

So I mounted the Trimble receiver onto a spare piece of PCB stock, added a switching voltage regulator, and connected the Trimble receiver's 1 PPS signal into my GPSDO in lieu of the Quectel L76 1 PPS signal.

The Trimble Resolution T was placed into "O-D Clk Mode" (i.e. Over-determined-Clock Mode) and set to "Stationary" Dynamics.

Here's how the receiver's 1 PPS jitter translates to phase-delta jitter in the model.  This was measured exactly the same way that the Quectel phase-delta jitter was measured (described in this post, earlier).

Note that the only thing I've changed in the model between this capture and the earlier capture from the Quectel L76 is the name of the Excel file.  In this case, the file's name is Trimble_PD_181227.xlsx.

There is quite a bit of peak-to-peak variation, but if we run it through the first IIR filter, the VCO response smooths out quite nicely:

Compare the above response to the earlier VCO response, using the same IIR filter, but with jitter from the Quectel L76 GPS receiver:

A significant difference!

This difference implies to me that the Trimble Resolution T receiver is probably a much choice of receiver to use for a phase-locked loop GPSDO application.

Note, thought, that with my own frequency-locked loop GPSDO, the Quectel L76 jitter performance has not been an issue.

By the way, the Trimble's jitter performance might not always be milk and honey.  Here's a capture I made the day after I made the plots, above:

That's a jump in my phase-detector's output of about 250 ADC counts.  And given that each count represents a time interval of 0.5 ns, this is a jump of about 125 ns!

Here is a zoomed-in look...

Of course, this phase jump could have been produced by my HP 10811-60111 Oscillator Module (which is my GPSDO's VCO), even with its EFC held at a constant value.  Unfortunately, it's not possible to determine which of the Phase Detector's two signals (1 Hz or 2.5 MHz) had the perturbation on it.


5.  Simulating a 5 MHz versus a 10 MHz VCO:

The design in Shera's original article uses a 5 MHz VCO that is divided down by 16 to create the 312.5 KHz signal fed to the phase detector.

The Simulink model described above, in this post, uses a 10 MHz.  To model a 5 MHz VCO instead, simply change Kv from 0.075 to 0.0375 (to maintain Shera's sensitivity 'S' factor of 7.5E-9) and change the divide factor N from 32 to 16.

Here's how the Simulink VCO response looks for 10 MHz versus 5 MHz:

The amplitude of the 5 MHz VCO response is half that of the 10 MHz VCO's response, but the settling time is identical for both.

The DAC input (and thus the DAC output) for both the 10 MHz simulation and the 5 MHz simulation are identical.


6.  Brook Shera's Digital IIR Filter:

From earlier in this post, Shera's equation for his IIR filter is:

o(n) = o(n-1) + i(n)*(1/F1 + 1/F2) + i(n-1)*(1/F1 - 1/F2)

In general, an IIR digital filter's output sequence, as a function of its input sequence, can be written as:

Comparing this equation with Shera's and recognizing that his o(n) is the same as y(n) and that his i(n) is the same as x(n), the terms for the second equation become:

• M = 1
• N = 1
• a₀ = 1/F1 + 1/F2
• a₁ = 1/F1 - 1/F2
• b₁ = -1

The general form of the IIR equation, above, has the following transfer function in the z domain:

In Shera's case, this transfer function becomes:

(For more on digital filter equations and transfer functions, see: Harry Y-F. Lam, Analog and Digital Filters: Design and Realization, Englewood Cliffs, N.J.: Prentice-Hall, Inc., 1979)


Simulink Model Repository:

My Simulink models (and other information) can be found at the following GitHub repository:
https://github.com/k6jca/Brooks-Shera-GPSDO


My  GPSDO posts:

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

(This post) 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 Phase-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, who is an expert in MATLAB and Simulink and whose advice and recommendations helped me break through a number of mental log-jams blocking my progress.


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.

And so I should add -- this 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.