Broadband Fringe-Fitting Algorithm

This section describes the core single-baseline, single-polarization-product broadband (bandwidth synthesis) fringe-fitting algorithm as implemented in HOPS4 (ported from HOPS3). This algorithm performs a three-dimensional search over single-band delay (SBD), multi-band delay (MBD), and delay rate (DR) to determine the residual geometric delay parameters of a VLBI observation after initial correlation. The search employs a coarse FFT-accelerated grid search followed by a fine sub-bin interpolation step to achieve sub-sample precision. The program which implements this algorithm (and other associated calibration operations) is called fourfit, or specifically in HOPS4 (fourfit4).

Introduction

In Very Long Baseline Interferometry (VLBI) each radio telescope (a station) in the array records time-tagged samples of a noise-dominated voltage signal collected from a common astronomical source. The resulting signals from each station are then cross-correlated with an appropriately tuned (a priori) delay model in order to compute the complex visibilities of the source over each baseline (station pair). The first station in a baseline pair, A, is referred to as the reference station, while the second station, B, is know as the remote station [1]. The resulting visibility data contains residual geometric phase terms due to imperfect knowledge of the source position, station coordinates, clock offsets, and atmosphere. Fringe fitting is the process of searching for the set of residual corrections to the delay model parameters that maximize the fringe amplitude which is produced from the coherent combination of the visibility data across all spectral channels and time intervals of an observation.

In HOPS4, fourfit4 implements a broadband (or bandwidth-synthesis) fringe-fitting algorithm that operates on the full multi-channel visibility data. This algorithm was first described in [Rogers1970]. This broadband fringe fitting process simultaneously resolves two main residual corrections to the delay model for a particular baseline. These are the multi-band delay (MBD, otherwise known as the group-delay), and the delay rate (DR, the time derivative of the group-delay).

In HOPS4 this algorithm is implemented across several source files:

  • Fringe/MHO_BasicFringeFitter.cc – main orchestration and pipeline control

  • Fringe/MHO_IonosphericFringeFitter.cc – extended fitter with dTEC correction/search

  • Calibration/MHO_SingleSidebandNormFX.cc – normalization and frequency-to-lag (SBD) transform (for USB or LSB only channels)

  • Calibration/MHO_MixedSidebandNormFX.cc – normalization and frequency-to-lag (SBD) transform (for mixed LSB/USB/DSB channels)

  • Calibration/MHO_MBDelaySearch.cc – coarse three-dimensional grid search over (SBD, MBD, DR)

  • Calibration/MHO_InterpolateFringePeakOptimized.cc – fine sub-bin fringe peak interpolation

  • Fringe/MHO_ComputePlotData.cc – post-solution diagnostics and visualization quantities

Input Data Model

In this document we will restrict the discussion to the visibility data, but for a complete specification of the data structures used by the fringe-fitting algorithm and HOPS4 see: Data Object Specification.

Visibility Tensor

The fundamental input is a four-dimensional tensor of complex visibilities, \(\mathbf{V}\), associated with each baseline:

\[V[p,\,c,\,a,\,s] \in \mathbb{C}\]

where the axes are:

  • p POLPROD_AXIS: polarization product index (XX, XY, RR, LR, etc.).

  • c CHANNEL_AXIS: channel index, this is the coarse frequency subdivision of the data. Each channel carries metadata: sky frequency \(\nu_c\) (MHz), net sideband (USB, LSB, or DSB), and bandwidth \(B_c\) (MHz).

  • a TIME_AXIS (or accumulation period (AP)): time/AP index, each AP spans a fixed integration time \(\Delta t\) (typically 1 s). The units of this axis are time (sec).

  • s FREQ_AXIS: frequency spectral-bin (or spectral point) within each channel (typically there are 64-512 bins, depending on the setup at the time of correlation). Units are MHz.

Note

The fringe-fitting algorithm operates on only a single polarization product (or linear combination of polarization products) at one time, and after user selection, the index p is essentially fixed with p=p’=0.

Weight Tensor

Each visibility is accompanied by a non-negative real weights tensor, \(\mathbf{W}\):

\[W[p,\,c,\,a,\,s] \in \mathbb{R}_{\geq 0}\]

The weights encode the data quality (e.g., zeroed for invalid/missing data) reported by the correlator. However, the values of the weights may be further modified by additional data operations (i.e. flagging) before fringe fitting.

Note

In fourfit4 the size along the 4th (FREQ_AXIS), s, dimension is always set to 1, fourfit4 does not yet support intra-channel weights for individual spectral points.

In order to form weighted averages over the visibilities, the total weights for the selected polarization product is calculated as:

\[W_{\mathrm{tot}} = \sum_{c,a} W[p',c,a,0]\]

Step 1: Normalization and Frequency-to-Lag Transform

The first computational step in the fourfit4 fringe fitting algorithm, transforms the visibility data from the frequency domain into a single-band delay (SBD) domain (indexed by lag). This operator, called NormFX, is implemented in MHO_SingleSidebandNormFX.cc for single-sideband data and MHO_MixedSidebandNormFX for mixed (USB/LSB/DSB) data. This transformation is applied only along the intra-channel FREQ_AXIS (indexed by spectral-bin s) of the visibility tensor.

Ideally, for intra-channel visibility data without any non-linear phase corruption, this operator has the effect of placing all of the fringe ‘power’ of each channel into a single lag bin. This makes further manipulation of the data during the (MBD, DR) search much more efficient. However, it must be noted that all channels must be phase and delay corrected such that they are aligned and share the same lag (SBD) bin. Otherwise, channels for which power is distributed to a different SBD bin will be lost, and not contribute coherently to the total fringe amplitude.

Zero-Padding

Let \(N_s\) denote the number of frequency spectral-bins. Then before the frequency to lag transformation, a padding factor of \(P = 4\) is applied: the FREQ_AXIS is extended to \(P \cdot N_s\) elements by appending zeros at the end. Zero-padding provides additional interpolation points in the lag domain, improving the resolution of the SBD peak search, and reducing the power which may bleed into other bins.

Forward Discrete Fourier Transform

A one-dimensional forward DFT is applied along the (zero-padded) frequency axis for each channel slice (p, c, a) of the visibilities, \(\mathbf{V}\). For a each channel c with bandwidth, \(B_c\), at AP a, we apply:

\[\mathcal{S}[p,c,a,\ell] \;=\; \sum_{s=0}^{N_s-1} V[p,c,a,s] \; \exp\!\left(-\frac{2\pi i \, \ell \, s}{P \cdot N_s}\right) \qquad \ell = 0,\dots,P\cdot N_s-1\]

The result \(\mathcal{S}\) is the SBD array (short for “single-band delay”). Each lag \(\ell\) corresponds to a time delay which is common to all channels. The full search range (window) of the SBD is determined by the channel bandwidth and the number of spectral points, \(N_s\), chosen during correlation, and is given by \(\pm \frac{N_s}{2 B_c}\). The delay resolution of the SBD axis is determined by the channel bandwidth and the padding factor (\(P = 4\)), and is given by \(\delta_{\mathrm{SBD}} = \frac{1}{P B_c}\). Note that this stage of the search algorithm requires all channels to have the same bandwidth, or to have been padded out to the same bandwidth, so that the SBD axis can be aligned across all channels.

Cyclic Rotation

The FFT output places the zero-delay (DC) term at index \(\ell=0\). To center the lag array around zero delay, a cyclic shift by \(P\cdot N_s / 2\) positions is applied along the FREQ_AXIS:

\[\mathcal{S}'[\ell] = \mathcal{S}\bigl[(\ell + P N_s/2) \bmod (P N_s)\bigr]\]

After rotation, negative delays appear in the first half of the array and positive delays in the second half, with zero delay at the center. Note that this operation is essentially the same as implemented by the matlab/python function fftshift.

Lower-Sideband Conjugation

Channels flagged as lower sideband (LSB, net_sideband = "L") have an inverted frequency axis relative to upper sideband (USB) data. To correct for this, the SBD slice for each LSB channel is complex-conjugated:

\[\mathcal{S}'[p,\,c_{\mathrm{LSB}},\,a,\,\ell] \;\leftarrow\; \overline{\mathcal{S}'[p,\,c_{\mathrm{LSB}},\,a,\,\ell]}\]

This is equivalent to flipping the frequency axis before the FFT, ensuring that USB and LSB channels contribute coherently to the same delay space.

Normalization

Finally, the entire SBD array is divided by the original (unpadded) number of frequency spectral-bins:

\[\mathcal{S}'[p,c,a,\ell] \;\leftarrow\; \frac{\mathcal{S}'[p,c,a,\ell]}{N_s}\]

This normalization compensates for the FFT summation over \(N_s\) terms. Neither the library FFTW3 nor the native FFT implementation apply this normalization by-default, so it must be applied after the FFT.

Output Dimensions

The SBD array, \(\mathcal{S}'\), has dimensions [1, N_c, N_a, P*N_s] (one polarization product, \(N_c\) channels, \(N_a\) APs, and \(P\cdot N_s\) lags). The labels associated with the 4-th axis are now physical delay values in microseconds.

Step 3: Fine Peak Interpolation

The coarse search yields integer bin indices. The fine interpolation (file: MHO_InterpolateFringePeakOptimized.cc) refines the peak location to sub-bin precision using a two-step process: (1) direct evaluation of the fringe function on a 5x5x5 grid, followed by (2) iterative 5-point Lagrange interpolation.

The 5x5x5 Grid

The grid is centered on the coarse peak \((\ell^*, k^*, m^*)\) with offsets:

  • SBD: 5 points at \(\ell^* - 2,\;\dots,\; \ell^* + 2\) (full-bin steps \(\delta_{\mathrm{SBD}}\))

  • MBD: 5 points at \(m^* - 2,\;\dots,\; m^* + 2\) (half-bin steps \(0.5 \cdot \delta_{\mathrm{MBD}}\))

  • DR: 5 points at \(k^* - 2,\;\dots,\; k^* + 2\) (half-bin steps \(0.5 \cdot \delta_{\mathrm{DR}}\))

This requires 125 trial points to fill the interpolation grid.

Evaluating Each Grid Point

For each trial point \((\ell_{\mathrm{trial}},\, m_{\mathrm{trial}},\, k_{\mathrm{trial}})\), the algorithm computes the weighted, counter-rotated coherent sum over all channels and APs:

\[Z = \frac{1}{W_{\mathrm{tot}}} \sum_{c=0}^{N_c-1} \sum_{a=0}^{N_a-1} W[p',c,a,0] \; \mathcal{S}'[p',c,a,\ell_{\mathrm{trial}}] \; \Psi(c, a, \tau_{\mathrm{mbd}}, \dot{\tau}_{\mathrm{dr}})\]

The fringe rotation phasor \(\Psi\) removes the residual delay and delay-rate terms from the visibilities:

\[\Psi(c, a, \tau_{\mathrm{mbd}}, \dot{\tau}_{\mathrm{dr}}) = \exp\!\Bigl(-2\pi i \bigl[ \nu_c \, \dot{\tau}_{\mathrm{dr}} \, \Delta t_a + \tau_{\mathrm{mbd}} \, (\nu_c - \nu_{\mathrm{ref}}) \bigr]\Bigr)\]

where:

  • \(\nu_c\) is the sky frequency of channel c (MHz)

  • \(\dot{\tau}_{\mathrm{dr}}\) is the trial delay rate (\(\mu\mathrm{s}/\mathrm{s}\))

  • \(\Delta t_a = t_a + \Delta t/2 - t_{\mathrm{FRT}}\) is the time offset of AP a from the fourfit reference time (FRT), in seconds

  • \(\tau_{\mathrm{mbd}}\) is the trial multi-band delay (us)

  • \(\nu_{\mathrm{ref}}\) is the reference frequency (MHz)

Precomputation and Phasor Recurrence

To avoid computing expensive \(\exp(\cdot)\) in the inner loop, the algorithm precomputes:

  • Channel phasors \(\Phi_c^{\mathrm{MBD}}[m_{\mathrm{trial}}]\) for each MBD trial point and channel (the frequency-offset part of the phasor)

  • Step phasors \(\Phi_c^{\mathrm{step}}[k_{\mathrm{trial}}] = \exp(-2\pi i \, \nu_c \, \dot{\tau}_{\mathrm{dr}} \, \Delta t)\) for each DR trial point and channel

  • Initial phasors \(\Phi_c^{\mathrm{init}}[k_{\mathrm{trial}}] = \exp(-2\pi i \, \nu_c \, \dot{\tau}_{\mathrm{dr}} \, \Delta t_0)\)

The phasor for AP a is then advanced by recurrence:

\[\Psi_c(a+1) = \Psi_c(a) \cdot \Phi_c^{\mathrm{step}}\]

This replaces \(N_a\) complex exponential evaluations with \(N_a\) complex multiplications.

Five-point Lagrange Interpolation

After filling the 5x5x5 amplitude cube:

\[\mathcal{F}[i,j,k] = |Z_{i,j,k}| \qquad i=0..4,\; j=0..4,\; k=0..4\]

the algorithm performs an iterative refinement search (function max555), which:

  1. Initializes the search location to the cube center, and step size dx = 0.4, in each dimension.

  2. Defines an 11x11x11 sub-grid within the current bounds.

  3. For each of the 1331 sub-grid points, computes the interpolated value using five point (quartic) Lagrange interpolation [AS1965] (25.2.15).

  4. Locates the maximum value and its coordinates.

  5. Relocate the center to the maximum, then reduces dx by a factor of 5.

  6. Repeats until dx < 1e-4 in all dimensions.

The three-dimensional interpolated value of the fringe amplitude at \((x_0, x_1, x_2)\) is given by:

\[\mathcal{F}_{\mathrm{interp}}(x_0, x_1, x_2) = \sum_{i=0}^{4}\sum_{j=0}^{4}\sum_{k=0}^{4} a_i^{(0)} \, a_j^{(1)} \, a_k^{(2)} \, \mathcal{F}[i,j,k]\]

using the following definition for the one-dimensional Lagrange coefficients at fractional position p:

\[\begin{split}\begin{aligned} a_0 &= \frac{(p^2-1)\,p\,(p-2)}{24}, & a_1 &= -\frac{(p-1)\,p\,(p^2-4)}{6}, \\ a_2 &= \frac{(p^2-1)\,(p^2-4)}{4}, & a_3 &= -\frac{(p+1)\,p\,(p^2-4)}{6}, \\ a_4 &= \frac{(p^2-1)\,p\,(p+2)}{24} \end{aligned}\end{split}\]

The final value for the peak location in interpolation space is given by:

\[(\xi_0, \xi_1, \xi_2) = \underset{(x_0,x_1,x_2)}{\arg\max} \; \mathcal{F}_{\mathrm{interp}}(x_0, x_1, x_2)\]

Refined Peak Parameters

The final refined (but dimensionless) coordinates \((\xi_0, \xi_1, \xi_2)\) are converted to physical units as follows:

\[\begin{split}\tau_{\mathrm{SBD}} &= \mathrm{SBD}[\ell^*] + \xi_0 \cdot \delta_{\mathrm{SBD}} \\ \tau_{\mathrm{MBD}} &= \mathrm{MBD}[m^*] + \xi_1 \cdot 0.5 \cdot \delta_{\mathrm{MBD}} \\ \dot{\tau}_{\mathrm{DR}} &= \mathrm{DR}[k^*] + \xi_2 \cdot 0.5 \cdot \delta_{\mathrm{DR}} \\ f_{\mathrm{rate}} &= \dot{\tau}_{\mathrm{DR}} \cdot \nu_{\mathrm{ref}} \\ A_{\mathrm{fringe}} &= \mathcal{F}_{\max} \quad \text{(interpolated peak amplitude)}\end{split}\]

Step 4: Compute Output

Once the location of the fringe peak has been determined. A collection of summary quantities and plots are calculated in order that the user may evaluate the quality of the fit. See the section Solution Quantities and Diagnostics for details.

Other Implementation Details

Ionospheric Fringe Fitting

For VGOS geodetic VLBI removing the effect of the ionosphere is an essential step during fringe-fitting. This modified algorithm is implemented in MHO_IonosphericFringeFitter. See the section Ionospheric Fringe Fitting: Solving for Differential TEC for details.

Parallel Backends

The coarse MBD search supports three parallel backends selected at compile time:

  • Sequential (MHO_MBDelaySearch) – single-threaded (default)

  • OpenMP (MHO_MBDelaySearchOpenMP) – shared-memory parallelism over the SBD lag loop

  • CUDA (MHO_MBDelaySearchCUDA) – GPU acceleration

FFT Backend

When FFTW3 is available, the algorithm uses the FFTW library for all DFT transformations. Otherwise, an internal FFT implementation (which may be slower) (MHO_MultidimensionalFastFourierTransform) is used. If CUDA is enabled, CUFFT is used, but only during the coarse MBD search.

Mixed-Sideband Handling

When the data contains a mixture of USB and LSB channels (or DSB channel pairs), the algorithm uses MHO_MixedSidebandNormFX instead of MHO_SingleSidebandNormFX. The mixed-sideband version applies frequency-axis flipping and conjugation on a per-channel basis before the FFT, ensuring that all channels contribute coherently to the same delay space.

Single-Channel Degeneracy

When only one active channel is present, the MBD is set equal to the SBD (by definition, a single channel cannot distinguish between SBD and MBD).

Caching and Iteration

The fringe fitter supports iterative refinement loops controlled by user-specified prefit and postfit operators. When operators from both categories are present, the original visibility and weight data are cached and can be refreshed between iterations, allowing convergence-based outer loops (e.g., iterative flagging or passband estimation).