1095 REAL–TIME PROCESSING OF DOWNLINK PHASE MEASUREMENTS

Larry R. D’Addario

National Radio Astronomy Observatory* 2015 Ivy Road, Charlottesville, Virginia 22903, U.S.A.

Introduction

The two-way time transfer link between earth station and satellite will be of the form shown in Figure 1. The downlink signal will be converted to a convenient I.F. and then its phase will be measured with respect to a programmable reference. Some details of the analog signal processing leading up to the phase detector were discussed in OVLBI-ES Memo 22. In our implementation, the phase detector is an analog device consisting of two mixers operating in phase quadrature. Outputs of the two mixers will be sampled and digitized at a rapid rate and the resulting data will be processed in real time by a digital signal processor. The processed data will be presented at a much slower rate to the main station control computer, which could perform additional real-time processing if necessary, but which is responsible mainly for writing the data to a log file. Off-line (non-real-time) software will convert the log file to the forms needed by correlators and orbit determination centers.

Fig. 1: Conceptual block diagram of the two-way time transfer.

The present memo is concerned with the design of the required real-time digital signal processing.

Sampling Rates

As mentioned in OVLBI-ES Memo 22, the highest feasible input sampling rate will be supported subject to a minimum requirement of 1 kHz. The calculations below will consider sampling rate of 10 kHz. At the DSP output (to the station computer), we assume a rate of 10 measurements per second. This is expected to be adequate for representing the fastest fluctuations in the signal phase, which are expected to be due to fluctuations in the tropospheric delay.

Typeset by AMS-TEX

*The National Radio Astronomy Observatory (NRAO) is operated by Associated Universities, Inc., under a cooperative  agreement  with  the  National  Science  Foundation.

Processing of Individual Samples

Certain calculations must be done for each individual sample:

1. Accumulate data on phase detector errors

2. Correct for phase detector errors

3. Convert measurement from rectangular (x, y) to polar (a, ϕ) form

4. Resolve cycle ambiguity (i.e., determine whether a cycle has been completed since the preceding sample, and if so increment or decrement the integer part of the phase)

The result of these calculations will be a single number equal to the measured phase. It needs to have a resolution of better than 0.1 radian (.016 cycle, 6 bits after binary point) to avoid significant phase noise being added, and a range sufficient to cover the maximum variation in orbital position error during a tracking pass. The latter is predicted to be < ±1 km, but larger errors might occur in the early phases of the missions. At a wavelength of 2 cm (VSOP), ±1 km gives a two-way phase error of ±10^{5} cycles, requiring 18 bits for the integer part of the phase. Thus, a 24-bit fixed point representation would be adequate but additional bits would be desirable.

The analog phase detector operation may be modeled as

\begin{aligned}x & = g_{x}a cos(ϕ) + x_{0} \\ y & = g_{y}a sin(ϕ + ϵ) + y_{0} \end{aligned}  (1)

where x, y are the phase detector outputs; x_{0}, y_{0} are d.c. offset errors, g_{x}, g_{y} are scale factors; and ϵ is the error in phase quadrature between the x and y measurements. The signal being measured has amplitude a and phase ϕ. Equations (1) are easily solved for a and ϕ:

x̂ = (x − x_{0})/g_{x}
ŷ = [(y − y_{0})/g_{y} − x̂ sin ϵ]/ cos ϵ
a = \sqrt{x̂^{2} + ŷ^{2}}
ϕ = tan^{−1}(ŷ, x̂)

(2)

where the arctangent function is the four-quadrant version (2π radian range, ATAN2 in Fortran). The parameters x_{0}, y_{0}, g_{x}, g_{y}, and ϵ can be determined by separate calibration measurements using a known test signal, but they can also be derived from a series of measurements on the unknown signal provided that they are relatively constant and that the signal amplitude a is also relatively constant. The process is simplified if ϵ = 0 to sufficient accuracy. In our implementation, the quadrature reference signals to the phase detector have separately programmable phase shifts which should allow us to make ϵ = 0; tests will determine whether this is accurately achieved over a range of signal frequencies and over temperature and time. To allow for determining the calibration parameters, the DSP will accumulate the sum of x, y, x^{2}, y^{2}, and xy over a programmable interval. These will be available to the station computer.

Resolution of the cycle ambiguity can be accomplished as follows. For each sample, let −π < ϕ ≤ π, and let the previous sample value be ϕ_{old}. Then calculate δϕ = ϕ − ϕ_{old}; if δϕ > π then decrement the integer part of the phase; if δϕ < −π then increment the integer part of the phase; otherwise leave the integer part unchanged. Note that in the presence of noise this algorithm can sometimes make the wrong decision. In that case, an unrecoverable cycle slip occurs. The probability of such a cycle slip must be kept extremely low, and the required signal-to-noise ratio depends on the ratio of the signal frequency to the sampling rate; the appropriate theory is developed in Appendix A. We proceed by assuming that the SNR is adequate.

The above algorithm becomes particularly simple when ϕ is represented in units of cycles using fixed point, twos-complement binary with the binary point just after the sign bit. Then the integer part is changed if and only if the δϕ subtraction causes an overflow; and the decision to increment or decrement depends only on the sign of δϕ.

Filtering

The above processing may result in 1000 to 10 000 phase measurements per second, each being a 24-bit (or larger) number. Once the phase ambiguity has been resolved, it is believed that the intrinsic bandwidth in this time series is no more than a few Hz. This is true in spite of the fact that over an interval of a 1 sec the measured phase may change by hundreds of cycles; this large time derivative can result from a component (viewed in the frequency domain) whose amplitude is around 10^{5} cycles but whose frequency is the order of 10^{−4} Hz. Such components are expected to result from position errors that vary on time scales related to the orbital period. If the signal is passed through a linear lowpass filter whose bandwidth is several Hz, components at such low frequencies will be negligibly affected, despite their large amplitudes, provided only that numerical errors (roundoffs and overflows) are properly controlled.

We therefore plan to process the individual high-rate samples, after cycle resolution, by passing them through a digital, finite-impulse-response lowpass filter of bandwidth 5 Hz and then re-sampling at 10 Hz. An FIR filter with symmetrical coefficients is always a constantdelay filter, which is a desirable feature in our application. Filtering and resampling are efficiently implemented together, forming what is known as a decimation filter.

Suppose that the filter is required to have less than 1.0% ripple in a passband of 0 to 4.0 Hz, and at least 40 dB of rejection in a stopband beginning at 5.0 Hz. Then to decimate from a rate of 1 kHz to 10 Hz (decimation factor M = 100) requires a filter length of about N = 1800 coefficients. With similar specifications, decimation from 10 kHz to 10 Hz (M = 1000) would require N ≈ 18000. Such long filters are impractical to implement; they require too much memory and computing time, and would accumulate excessive numerical errors. However, there is a way to implement the filter in several stages so as to avoid these problems [1]. Consider the following three filters in cascade. For each one, the passband extends from 0 to 4.5 Hz with ripple less than 0.5%, and the stopband rejection is 40 dB.

I. Decimate 10 kHz to 500 Hz, stopband starts at 495 Hz; needs N_{1} = 49 coefficients. II. Decimate 500 Hz to 50 Hz, stopband starts at 45 Hz; N_{2} = 39. III. Decimate 50 Hz to 10 Hz, stopband starts at 5 Hz; N_{3} = 115.

The output of each stage is the input of the next. The overall result is a decimation from 10 kHz to 10 Hz with the desired response, and a total of N_{1} + N_{2} + N_{3} = 203 coefficients is needed (vs. 18000 with a single stage). We will use this approach in our implementation.

Signal Processors

We have so far been experimenting with the Analog Devices ADSP-2100 series of digital signal processor chips; in particular, we are using the ADSP-2105 in our laboratory breadboard. The overall architecture of this processor is shown in Figure 2. This did not result from any rigorous selection process, but rather from the fact that a single-board computer using this processor was available at low cost along with some development software. The device appears to be adequate for our purposes, although the choice is not final.

Fig. 2: Block diagram of the ADSP-2105 processor (from [2]).

The ADSP-2105 is a fixed-point processor with 16-bit data words and 24-bit program words in separate memories. All instructions execute in exactly 100 nsec (1 cycle of the external clock), and calculations can be pipelined so that several things happen in the same cycle. For further details, consult the manufacturer’s literature [2,3]. Calculations of computing times given below are based on the use of this processor.

Implementation Details

. Hardware

A block diagram of the prototype hardware configuration is shown in Figure 3. Timing of the entire subsystem is controlled by the digital signal processor; its 10 MHz clock is taken from the station timing reference (H maser), and it receives 1 Hz interrupts synchronized to UTC seconds. The DSP contains a programmable timer that enables any frequency between these to be synthesized. In particular, the sampling frequency is generated and used to trigger conversions in two analog to digital converters connected to the (filtered) phase detector outputs.

The DSP also controls the two direct digital synthesizers. Present plans are to update both of them once per downlink phase sample; although this is more often than necessary, the CPU cycle requirements are modest. Further discussion of the DDS programming will be given in a separate report.

Communication with the station computer is via the VLBA Monitor/Control Bus (like all of the station hardware). A shared memory of 16-bit words is addressable by both the MCB and the DSP, and is used to pass control parameters to the DSP (mainly the DDS frequency parameters) and to pass the downlink phase measurement results to the station computer. Diagnostic information and calibration results are also passed. The DSP operates most of the time with instructions and data stored in on-chip memory, so that access to the shared memory by the MCB does not affect execution of DSP instructions. But if a DSP instruction requires access to external memory, including the shared memory and the ADCs (which are memory mapped), it must wait during any MCB access.

To ensure precise timing of critical events like the phase detector sampling and the DDS updating, MCB access can be inhibited under software control. With this precaution,

Fig. 3: Digital subsystem block diagram.

and considering that all DSP instructions execute in exactly one cycle and that the code is deterministic, the sampling times and DDS update times are predictable with a resolution of one DSP clock cycle (100 nsec) and a precision limited only by the clock jitter (probably a few nsec).

The DSP excecutable code is normally stored in an EPROM and is automatically loaded into the on-chip program memory during the power-up sequence. Under software control, the processor can be re-loaded with different code from the EPROM. Our configuration also makes it possible to download code from the station computer over the MCB.

. Firmware

An outline of the required code is given in Table 1. Calculations and simulations show that there should be little difficulty in completing the required computations in the available time, even for the 10 kHz sampling rate. Here are some of the timing estimates:

Operation  CPU Cycles Samples/Oper.
Trigger next sample; read ADCs  6  1
Accumulate calibration data  28  1
Correct for phase detector errors  11  1
Arctangent (maximum)  95  1
Resolve cycle wrap (maximum)  10  1
Update uplink and downlink synthesizers  67  1
Misc. flow control and test instructions  15  1

Filter stage 1

Copy data to buffer  4M_{1} + 7  M_{1}
Compute output  3N_{1} + 19  M_{1}
Filter stage 2—compute output  3N_{2} + 19  M_{1}M_{2}
Filter stage 3—compute output  3N_{3} + 19  M_{1}M_{2}M_{3}

The M_{i} are the decimation factors and the N_{i} are the lengths of the filter stages. The filter computation cycle counts are based on double-precision calculations (32-bit numbers) using code supplied by Analog Devices [3]. The other times are based on counting instructions in actual code (preliminary, untested version). These numbers are believed to be quite accurate.

As discussed earlier, we choose M_{1} = 20, M_{2} = 10, and M_{3} = 5 to achieve a total decimation of 1000 from 10 kHz to 10 Hz. Because of the memory addressing facilities built

6 Larry R. D’Addario

Table 1: Outline of DSP Firmware into the ADSP-2105, it is efficient if the filter lengths are of the form 2^{n} − 1; so we choose N_{1} = 63, N_{2} = 31, and N_{3} = 127 rather than 49, 39, and 115, respectively, which would just meet the specifications.

At the 10 kHz rate, there are 1000 CPU cycles available per sample (100 nsec/cycle). Using the values of M_{i} and N_{i} just selected, we find that at most 232 cycles are needed for the code that executes every sample, plus 255/20 cycles per sample for the first stage filtering, and 112/200 and 112/1000 cycles per sample for the respective other stages. Thus a total of 246 cycles per sample will be used, or less than 25% of the available cycles.

Present Status

All of the digital hardware shown in Figure 3 exists in breadboard form and has been tested successfully. Some experimental code for the DSP has been written and tested, including nearly all of the processing that must be done for each sample. The three-stage FIR filter has not yet been coded; a single-stage version has been written but not yet tested.

REFERENCES

[1] R. E. Crochiere and L. R. Rabiner, Multirate Digital Signal Processing. Prentice-Hall, 1983.

[2] Analog Devices Corp., “DSP microcomputer: ADSP-2105.” Data sheet, October 1990 (P.O. Box 9106, Norwood, MA 02062; 617/329-4700).

[3] Analog Devices Corp., A. Mar, ed., Digital Signal Processing Applications Using the ADSP-2100 Family. Prentice-Hall, 1990.

Appendix A: Required Signal-to-Noise Ratio

Let ϕ = (ϕ_{t} + e) mod 2π be the result of the most recent phase measurement, where ϕ_{t} is the true value of the phase and e is the measurement error (due to system noise). Similarly, let ϕ′ = (ϕ′_{t} + e′) mod 2π be the result of the previous phase measurement. Then the change in measured phase is δϕ = ϕ − ϕ′ and this quantity will be used to determine whether a cycle has been completed according to the algorithm

if δϕ > π then i := i − 1;

if δϕ < −π then i := i + 1 (A1)

where i is the cycle count. The total measured phase is then

θ = 2πi + ϕ.

We say that a cycle slip occurs if θ differs by more than π from the value it would have had if there were no noise (e = e′ = 0). The algorithm is based on the assumption that the true value of θ will never change by more than π from one sample to the next; this will be the case for sampling at or above the Nyquist rate in the absence of noise. A cycle slip occurs when the added noise causes the measured data to violate this rule. Thus, a cycle slip occurs whenever

|(θ − θ′) + (e − e′)| > π  (A2)

which implies

|θ − θ′| + |e − e′| > π  (A3)

and

|e − e′| > π − |θ − θ′|  (A4)

To limit the expected rate of cycle slips to 1/T, we require

Pr{|E| > π − |δθ|} < 1/(f_{s}T)  (A5)

where E = e − e′, δθ = θ − θ′, and f_{s} is the sampling rate. This assumes that all errors are independent and identically distributed and that δθ is nearly constant (i.e., the signal is quasi-sinusoidal). (For example, at f_{s} = 1 kHz and T = 1 hour, the probability must be < 2.78 × 10^{−7}.)

If the signal whose phase is being measured is

a sin(2πft) + n(t)

where n(t) is a zero-mean, stationary random process with ⟨n^{2}⟩ = σ^{2}, then any one measurement of the phase will have an error whose distribution is given approximately (for small errors) by

Pr{|e| > ϵ} = 0.5 [2 − erfc(a/σ)] erfc((a/σ) tan ϵ)  (A6)

where erfc(·) is the complementary error function. (The proof of this is omitted here for brevity, but can be supplied by the author upon request.) If the noise on successive samples is independent, then we can substitute E for e on the LHS of (A6) by substituting \sqrt{2}ϵ for ϵ on the RHS. Since δθ = 2πf/f_{s}, we have ϵ = π(1 − 2f/f_{s}). With these substitutions, and assuming that a/σ > 2, (A6) becomes

Pr\{cycle slip\} ≈ erfc[(a/\sqrt{2}σ) tan π(1 − 2f/f_{s})]  (A7)

and this can be evaluated for cases of interest. A few results follow.

\begin{aligned}(a/σ)^{2} = 10, f/f_{s} = 0.25 & ⇒ p = 1.34 × 10^{−23}. \\ (a/σ)^{2} = 10, f/f_{s} = 0.4 & ⇒ p = 0.22 \\ (a/σ)^{2} = 30, f/f_{s} = 0.4 & ⇒ p = 6.5 × 10^{−5} \\ (a/σ)^{2} = 100, f/f_{s} = 0.4 & ⇒ p = 3.3 × 10^{−13} \end{aligned}  (A8)

From this it can be seen that for signals whose frequency is less than half the Nyquist rate (f/f_{s} = 0.25), low SNRs (< 10 dB) are adequate; but as the Nyquist rate is approached things get rapidly worse, so that at 80% of Nyquist (f/f_{s} = 0.4) the SNR needs to be nearly 20 dB to keep the rate of cycle slipping to less than one per hour.