IWAVE—An adaptive filter approach to phase lock and the dynamic characterization of pseudo-harmonic waves

The co-inventor of the maser,¹ Arthur Schawlow, is said to have advised his students “Never measure anything but frequency!” There exist a great variety of techniques for measuring the frequencies of pseudo-harmonic waves, broadly addressing two different classes of problems.

In the first class, the average wave characteristics are estimated across a measurement interval under the implicit assumption that the wave properties are essentially static or when changes in these properties over the measurement interval are not of interest. This problem is addressed by a wide variety of methods, including Welch’s method using discrete Fourier transforms (DFTs),^2,3 Pisarenko’s method,⁴ MUSIC,⁵ and ESPRIT.⁶ 

In the second class, the oscillator is constantly evolving, and we seek a time-evolving best estimate of oscillator parameters. This problem is most often solved using phase locked loops (PLLs)⁷ or their hardware realization, lock-in amplifiers.⁸ The myriad applications of PLLs in science and engineering include, for example, a recent proposal for digital PLLs as an alternative technology for the readout of the photodiode clusters used to stabilize the alignment of mirrors in the future gravitational wave detectors.⁹ 

In searches for almost continuous wave (CW) signals in gravitational wave detectors, the evolving oscillation of continuous wave (CW) signals is at a very low signal-to-noise ratio (SNR), and so, the frequency evolution of such signals in ground based gravitational wave interferometers cannot be inferred from the raw data as is necessary for successful tracking with a conventional PLL. Instead, coherent matched filtering techniques, such as the $F$-statistic used in CW searches,^10,11 are able to achieve higher sensitivity to these weak signals at the price of substantially greater computational burden. Stack-slide-based semi-coherent algorithms expedite the computation to some extent at the cost of sensitivity by summing the signal power in multiple coherent segments after sliding the segments in the frequency domain to account for the signal phase evolution.^12–15 More efficient semi-coherent methods, e.g., signal tracking algorithms based on hidden Markov models,^16–19 have been developed to tackle the computational challenge as well as to allow for some uncertainties in the signal evolution model.^20–22 

Several other techniques beyond the conventional PLL have been developed for a range of applications outside the field of gravitational wave research although these algorithms are not adapted for the detection of the very weak CW signals expected from gravitational wave sources. One example is the second order generalized integrator (SOGI-PLL).²³ Developed for the problem of characterizing the characteristics of AC voltages in power lines, SOGI was one of the first methods to successfully address the problem of efficient generation of a copy of an input sinusoid that is out of phase (the so-called quadrature or Q phase) with the input signal. Other more recent papers, for example,^24–27 have developed the SOGI algorithm for a range of practical applications. A second example is the enhanced PLL (EPLL),^28,29 which solves the problem of tracking the amplitude of a harmonic wave in addition to its frequency. Further PLL developments are well summarized in Ref. 30.

The IWAVE (Iterative Wave Action-angle Variable Estimator) technique described in this paper is a new type of PLL addressing the dynamic characterization of evolving pseudo-sinusoidal signals. Unlike a conventional PLL, the adaptive element is a filter rather than an oscillator or counter. IWAVE has certain advantages over existing PLLs. First, IWAVE produces a benign output when the PLL is unlocked. In the case where the output of the PLL is being used to control something, for example, some parameter of a gravitational wave detector in a closed loop feedback system, then, in colloquial terms, IWAVE does no harm when it is not working. IWAVE naturally tracks the amplitude as well as frequency using a single feedback loop, unlike EPLL, which requires two loops internally for the amplitude and the phase. IWAVE is initialized using a small set of free parameters corresponding directly to physical oscillator properties—just an initial frequency and a single time constant; other PLL algorithms typically contain many control parameters that do not have a clear physical meaning. IWAVE also has the ability to characterize, simultaneously, multiple oscillations having almost-degenerate frequencies using the cross-subtraction method described in Sec. III C. This last advantage has led to detailed analysis of almost-frequency-degenerate violin modes of fused silica suspension wires in advanced LIGO.³¹ Comparisons of IWAVE with the SOGI and EPLL algorithms are given in the  Appendix.

As we discuss in Sec. III B, IWAVE as currently implemented is not sensitive to signals where the ratio of the amplitude of the target wave to the root mean square (rms) of the noise background is less than about 0.3. When applied to broadband gravitational wave data, therefore, IWAVE cannot be expected to detect gravitational wave CW signals at the strengths anticipated for sources described in the literature. However, preprocessing steps to divide the data into narrower frequency bands, thereby significantly reducing the rms of the noise, may yield promising search methods. Studies of these ideas are under investigation and will form the subject of future papers. The potential advantage is that by using the data itself to track the evolving frequency of the oscillation, IWAVE may be significantly less reliant on banks of templates for wave evolution and hence may require significantly less computational resources than existing CW search methods. For now, however, IWAVE represents a simple, well-characterized, and useful technique for analyzing quite complex spaces of evolving oscillators at relatively low computational burden, which is already being applied to studies of important background oscillations in gravitational wave data. This, then, is a method paper describing how IWAVE works and how it performs and giving some examples of that performance on gravitational wave data.

The IWAVE algorithm has many potential applications beyond gravitational wave science. In the control of brushless electric motors, IWAVE could out-perform standard vector control methods at low rotation rates where the back-emf in the motor windings is weak.³² In radio communications, IWAVE might be used to separate closely spaced channels with frequency evolution and multipath splitting.³³ In heart magnetometry, IWAVE could be used to estimate Fourier coefficients of the cardio-magnetic signal using an electrocardiograph signal to modify the phase evolution as the heartbeat rate evolves.³⁴ In physics, IWAVE could be applied to atomic force microscopy³⁵ or as part of an optical squeezing scheme for the readout of interferometers.³⁶ The authors look forward with anticipation to seeing what other applications we have not noticed.

This paper consists of a description of the IWAVE method in Sec. II, a discussion of the limits of applicability of IWAVE in Sec. III, and an overview of certain other PLL methods in the  Appendix. Space constraints have led us to leave out many mathematical steps; a full treatment of the mathematics can be found at Ref. 37. A software library implementing IWAVE in C with wrappers into MATLAB and Python/NumPY is available on a public git repository here.³⁸ 

Before writing down the IWAVE core algorithm, we consider as a starting point the $Z$-transform³⁹ of a regularly sampled time series, x_p,

where p increases with time and Ω = w − iΔ has the real part w the reciprocal of one e-folding for the weighting of previous samples and the imaginary part Δ equal to the frequency of the $Z$-transform component in radians per sample.

$Zn(Ω)$ obeys an iteration equation,

An unit amplitude phasor input, x_n = e^inΔ, to Eq. (2) results in the output $Zn=einΔ/(1−e−w)$⁠, which has zero phase shift with respect to the input and a larger amplitude dependent on w. Scaling the input by a factor of 1 − e^−w leads to an iteration algorithm, which passes phasors at frequency Δ with unit gain and zero phase shift,

This iteration algorithm is the core of IWAVE. As we shall see, it responds resonantly at frequency Δ. In the language of signal processing, Eq. (3) is an infinite impulse response (IIR) filter because it generates its nth output using the current input x_n, the previous output y_n−1, and a two-input, two-output, multi-input, multi-output (MIMO) filter; y_n and x_n are, in general, complex variables. We will also need the real representation of the transfer function for IWAVE derived from Eq. (3) by writing $xn=xnR+ixnI$⁠, $yn=ynR+iynI$⁠, and y_n−1 = z⁻¹y_n, where z⁻¹ is the sample delay operator. In these terms, Eq. (3) can be re-written as

where the elements of the transfer function matrix are

In this form, the transfer function is seen to be second order in sample delay. We determine the response of the core algorithm to an input consisting of a phasor of arbitrary frequency Θ radians per sample by substituting x_n = e^inΘ into Eq. (3). The response is also a phasor, $yn=AeinΘ+Φ$⁠, where $AΘ$ and $ΦΘ$ are the frequency dependent magnitude and phase lag of the output phasor with respect to the input given by

Thus, phasors of arbitrary frequency are eigenfunctions of the core algorithm with eigenvalues Ae^iΦ. The resonant character of the eigenvalues at frequency Δ can be seen in the denominator of Eq. (5), which is identical to the resonant denominator in the SOGI filter.²³ Figure 1 shows A and Φ as a function of Θ for Δ = 1.257 and four values of w: 1, 0.1, 0.01, and 0.001.

Starting from Eq. (6) and writing δ = (Δ − Θ) ≪ 1 and w ≪ 1 so that we are in a limit where the frequency is in the vicinity of a narrow resonance, the magnitude of the filter output can be approximated as

so that the peak is approximately Lorentzian in shape with full width at half maximum (FWHM) of 2w radians per sample. Using the sampling rate, f_s, in Hz, we give other properties of narrow resonances occurring when w ≪ 1 in Table I, where τ_s is the sampling period in seconds, τ is the response time in seconds, and Δ₀ is the resonant frequency in radians per sample.

We next discuss the usual case where the input data are a real oscillation at the IWAVE resonant frequency, $xn=cosnΔ$⁠. We decompose x_n into two phasors, x_n = x_f + x_b, where x_f = e^+inΔ/2 and x_b = e^−inΔ/2, each of which is an eigenfunction of the core algorithm. Substituting the frequencies ±Δ into Eq. (6) and rearranging, we obtain the response

where

By inspection, the locus of y_n is an ellipse in the complex plane having semi-major and semi-minor axes 1 + a and 1 − a. Figure 2 shows the result of driving the IWAVE core algorithm with an input x_n = cos nΔ for n ≥ 0. The output starts at the origin, spiraling outward toward a limiting ellipse in the steady state. Notice that the argument of the output y_n is always the same as the phase, nΔ, of the input sinusoid.

Input sinusoids at frequencies other than Δ also result in an elliptical steady state output although with different inclination angles and eccentricities and with smaller overall areas due to the falloff in the magnitude of the response for phasors having frequencies far from Δ.

A matrix transformation can be used to transform the elliptical locus into a circular one with the real part of this circular locus being a sinusoid in-phase with the input wave and the imaginary part having the same amplitude but lagging the input wave by 90°. We refer to these in-phase and out-of-phase components as the D and Q phases, respectively. This transformation can be expressed as a sequence of three elementary operations on the vector whose elements are the real and imaginary parts of y_n: a rotation through an angle of $ϕ2$ about the origin; shears parallel to the real and imaginary axes by factors of $21+a$ and $21−a$⁠, respectively; and finally a rotation through an angle of $−ϕ2$ about the origin. These three matrices can be combined into a single transformation,

The IWAVE core algorithm followed by this matrix transformation results in the output of both D phase and Q phase copies of the input drive. Thus, IWAVE is an example of what is referred to in signal processing parlance as an orthogonal state generator. Our convention follows other papers in the field in that the D (Q) phase output is in (out of) phase with the input at resonance. As we shall see, the Q phase quadrature can be used to generate an error signal to detect changes in frequency, allowing IWAVE to be used in place of a reference oscillator in a PLL. The sum in quadrature of the two phases, $An=Dn2+Qn2$⁠, is an estimate of the input signal amplitude. We will use the symbols in Fig. 3 to denote the application of IWAVE either to complex phasor or real sinusoidal inputs.

We next consider the transfer functions from a real sinusoidal drive, an example of which is shown in Fig. 4. Note that the dependence on frequency off-resonance is different for the D and Q outputs. The D transfer function rises linearly in frequency below the resonance and falls linearly in frequency above it, having a phase lead of 90° below the resonance and a phase lag of 90° above it. The Q transfer function has a flat frequency response below the resonance but falls as f⁻² above it and is in phase with the drive below the resonance but 180° out of phase with the drive above it.

This behavior is similar to that of a driven series RLC tank circuit, where the transfer functions from the input voltage to the voltage across the capacitor and resistor are similar to those between the input and the Q and D phase outputs, respectively. The relatively light suppression of low frequency off-resonance signals in the Q phase output can be important, particularly in cases where the quality factor Q_f of the circuit is set low by using a relatively large w coefficient. As in the resonant circuit, the ratio of the resonant to low frequency response is Q_f.

Changes in the amplitude of the incoming wave at the resonance result in a corresponding change in the quadrature sum, $An=Dn2+Qn2$⁠, but with a response time τ = τ_s/w leading to a single pole in the response to amplitude changes at s = −1/τ. Figure 5 shows the response at IWAVE’s resonance frequency to an input having a sinusoidally modulated amplitude for a variety of response times.

In order to phase lock with IWAVE, we need a measure of departures in the frequency from the frequency Δ of the harmonic wave at the input. We achieve this by exploiting the response time τ of the IWAVE algorithm. Consider a harmonic wave initially at frequency Δ. IWAVE yields both a D_n phase and a Q_n phase copy of this wave at its output. The product of the out-of-phase copy and the input wave, $A2⁡cosnΔsin(nΔ)=A2/2sin2nΔ$⁠, is a pure harmonic signal at frequency 2Δ. Now, consider an input signal where the wave develops anomalous phase and amplitude disturbances, δ and ɛ, respectively, so that x_n = A(1 + ɛ)cos(nΔ + δ). For elapsed times significantly less than τ following the onset of these disturbances, the outputs D_n = A cos nΔ and Q_n = A sin nΔ are unaffected by them. Physically, where the oscillator has an unchanged frequency and amplitude, its complex plane representation precesses in a circle about the origin, and the current IWAVE filter output therefore provides a predictor of the point where such a static oscillator will land the next sample. If the frequency or amplitude evolves, then the combinations E_n = (x_n − D_n)Q_n and $Fn=xnDn+Qn2−An2$ are estimates of the departure of the actual evolution of the complex coordinate of the oscillator state from this static model. These combinations can be written in a matrix form as follows:

This equation shows that E_n and F_n each consists of a static offset plus upper sidebands of frequency 2Δ, linear in the phase and amplitude offsets, respectively. The upper sidebands, however, appear as components of a rotating phasor in the space spanned by the E_n and F_n signals. We remove them using the complex carrier version of IWAVE, subtracting its outputs from its inputs. We have determined experimentally that using twice the w factor in the 2Δ IWAVE filter yields an acceptable error signal for the detection of phase departures.

Figure 6 is a schematic for the full IWAVE-based phase and amplitude detector. We have only considered here the case where the input data are real and we are tracking harmonic waves. The case where the input data are two-component rotating phasors is simpler as it can be shown that the combination $En′=xnRQn−xnIDn$ contains a pure DC signal in-phase offset with no upper sideband contamination. In addition, shown is the feedback path from the filtered error signal, δϕ_n, back to Δ through an integrator, discussed in the following.

The response of IWAVE to modulation of the frequency of the carrier is shown in Fig. 7. Knowing the analytic form of the frequency response, we can analyze the closed loop IWAVE PLL to determine an appropriate choice of feedback gain, G. A schematic for the feedback controller is shown in Fig. 8. Any difference between the incoming wave frequency and the IWAVE filter central frequency results in an accumulating phase shift in the homodyne detector. This accumulation of phase is represented by the factor of 2π/s, where s = 2πif with f being the signal’s frequency. The response of IWAVE to phase, confirmed by the measurements underlying Fig. 7, acts on the accumulated phase to produce the error signal. As shown in Fig. 8, the feedback path from the error signal to a correction in the central frequency of IWAVE consists of an adjustable gain, G, an addition of the gain boosted error signal to the previous value of the phase shift per sample, and a scale factor to convert from radians per sample to frequency in Hz. The closed loop gain is calculated by the usual consistency argument around the loop,

from which we obtain the closed loop transfer function,

This is the transfer function of a driven damped harmonic oscillator. We want a critically damped response since the control signal will be used to measure the oscillator frequency. For critical damping, we require two coincident real poles, which is achieved if $G=τs2/(4τ2)$⁠. The closed loop transfer function then takes the simpler form

The response to frequency or phase modulation is therefore flat below the knee frequency of 1/(4πτ) Hz, where due to the two poles it has rolled off to −6 dB and drops as 1/f² above that frequency. Larger values of G result in sharper turnover or a resonant peak, corresponding to the underdamped case.

The time constant, τ, determines the responsiveness of IWAVE to changes in wave frequency and amplitude as well as the noise bandwidth of the filter. Decreasing τ results in a faster response and better ability to stay locked on waves whose frequencies are changing but also increases the bandwidth of IWAVE to input noise, resulting in a noisier error signal. The optimal τ is small enough so that IWAVE stays locked when the frequency changes but not so small that excessive background noise is admitted by the filter. Section III A is on frequency tracking, and Sec. III B discusses locking in the presence of additive noise and the character of the error signal.

Consider a wave whose displacement at time t is $ht=A⁡cos2πftt$ so that the frequency of the wave is changing. The ability of IWAVE to track the evolving wave depends on the value of the response time parameter, τ. Figure 9 shows the results of running IWAVE on a swept sinusoidal wave starting at 20 Hz and increasing linearly in frequency at a variety of rates. At each sweep rate, a variety of values of τ were used. The sum of the squares of the deviation between the actual frequency and the frequency reconstructed by IWAVE over a time interval of 20 s, here referred to as χ², was calculated as a measure of the accuracy of frequency reconstruction.

The value of χ² is seen to be a function of τ with a pronounced minimum that is a function of the sweep rate. The χ² statistic also becomes noisy at both ends of the range of values of τ. Here, we explain the smooth descent and ascent in χ² on either side of the minimum and obtain a formula for the optimum value of τ. We also explain the onset of noise at either end of the range of τ.

Starting at the minimum, χ² rises with τ because of the response time of the closed loop transfer function of the servo from Eq. (14). This transfer function is that of two RC lowpass filters in series, each having an exponentially decaying impulse response of timescale 2τ. The response time of the entire transfer function is therefore the time duration of the autocorrelation of this exponentially decaying function. This autocorrelation initially rises linearly in time after the impulse before reaching a maximum at time 2τ and then decaying exponentially. The time between the impulse and the point where the exponential decay reaches 1/e of its maximum value is 6τ. This causes the frequency tracking of IWAVE to lag behind that of the input swept sine wave, leading to a frequency discrepancy at any given time of Δf₁ = 6τ(df/dt).

Below the minimum of χ²(τ), the rise in χ² with decreasing τ is explained by the frequency response of the core IWAVE algorithm and is best understood by the analogy between the IWAVE algorithm and the characteristics of a damped harmonic oscillator discussed in Sec. II B. At higher values of damping, the frequency response is maximal at a frequency f below the natural frequency of the undamped oscillator, f₀ by an amount given by the relationship $f0=f2+1/2π2τ2$⁠. This causes a systematic offset between the frequency, f, returned by IWAVE, and the frequency, f₀, of the input signal. The square of this frequency difference is an additional contribution to the χ² statistic.

By squaring and adding the frequency discrepancies arising from these two effects and minimizing with respect to τ, making the assumption that 2π²τ²(df/dt)² ≫ 1, we arrive at a value of τ that minimizes the χ² statistic for frequency tracking,

Figure 9 shows in bold red the values of τ_opt corresponding to the measured minimum τ and χ² across each of the χ² curves vs the measured χ² at that minimum along with in bold blue the value of τ_opt vs the value of χ² at τ_opt. There is good agreement between the theoretical τ_opt and the value determined from simulations within 29% for the largest value of df/dt studied and within 13% for the smallest one. At larger values of df/dt, τ_opt is smaller, so there is more broadband noise in the error signal, and a more sophisticated optimization on τ would lead to a larger optimal value than that predicted by Eq. (15).

We next discuss the breakdown of IWAVE at high values of τ, where the χ²(τ) curves become noisy, indicating loss of lock. The following argument leads to successful prediction of the value of τ where this breakdown occurs. Consider a wave whose displacement at time t is $h(t)=A⁡cos2πf(t)t$ so that the frequency of the wave is changing. Assume that IWAVE is locked at time t so that the IWAVE output is m(t) and is equal to h(t). At time t + τ, the wave displacement has evolved to $h(t+τ)=A⁡cos2π(f+df)(t+τ)$⁠. In the time interval [t, t + τ], the IWAVE output has not had time to respond to the frequency shift df and therefore takes the form $m(t+τ)=A⁡cos2πf(t+τ)$⁠. The phase shift between the wave and the IWAVE output at time t + τ is Δϕ ≃ 2πτ df. The error signal for IWAVE is approximated by the integral ∫h(t′)m(t′) dt′ over time interval [t, t + τ]. If the phase shift between the incoming wave and the IWAVE output over this time interval is greater than π, the error signal will undergo a sign change causing loss of lock. Therefore, a condition for IWAVE to remain locked is that Δϕ < π or τΔf < 1/2. Writing df = τ × df(t)/dt, we arrive at the upper bound that τ should obey

for IWAVE to remain locked. This line is drawn on the χ² curves in Fig. 9 in bold black, labeled U.L., vs the measured value of χ² at that value of τ in each simulation. The limit tracks the onset of lock loss as demonstrated by large excursions in χ² well across different trial values of df/dt.

The error signal for the IWAVE PLL is derived from the product of the input data stream and the Q phase output of the IWAVE filter minus the product DQ of the two IWAVE outputs. Where the wave frequency is static, this subtraction removes the upper sideband component at frequency 2f. When the frequency changes, this causes an additional transient upper sideband component, which is removed using a second IWAVE filter at frequency 2f having a response time half that of the primary IWAVE filter. Finally, we divide by the square of the wave amplitude because both the amplitude of the incoming wave and the amplitudes of both IWAVE outputs scale linearly with the wave amplitude. We need to divide this scale out; otherwise, the PLL loop gain will be dependent on the wave amplitude.

This is a form of homodyne detector. Because the error signal incorporates the unfiltered wideband input, the error signal incorporates the broadband noise of the incoming data. The distribution of the error signal therefore reflects the spectral characteristics of the input data over the full Nyquist band. If, for example, the incoming data includes a time domain transient with a broad spectral distribution, this transient will be reflected in the time history of the error signal from IWAVE. If the out of band noise is sufficiently large in amplitude, then IWAVE will lose lock.

We have determined experimentally that at the optimal choice of response time, τ_opt given in Eq. (15), then IWAVE will stay locked when the ratio of the wave peak amplitude to the root mean square noise amplitude exceeds 0.3. Future work to improve the performance of IWAVE at lower signal-to-noise ratios could involve, for example, prefiltering the input data to focus on a narrower frequency band about the frequency of interest for the waves under study. The noise content of the error signal, which leads to noise also in the estimate of the IWAVE frequency, is greater at smaller values of τ where the bandwidth of the IWAVE filter resonance is larger. At sufficiently small values of τ, the incursion of noise leads again to loss of lock. This can be seen in Fig. 9 for τ < 0.02. The optimal value of τ is affected by higher noise levels with larger values than that given in Eq. (15) becoming optimal.

The exact value of τ where lock loss occurs and the effects of noise on the optimal value of τ could be determined by a detailed stochastic differential equation analysis and is beyond the scope of this paper. However, a simple scaling argument can be made. Noise in the error signal leads to a random component being added to the phase shift per sample. This means that the reconstructed frequency, which also forms the servo control signal, contains a component that undergoes a random walk and hence grows with the square root of the number of samples. This means that you might expect the onset of loss of lock to occur at τ, which scales as one over the square of the signal-to-noise ratio, so that IWAVE works at τ above a lower limit that goes down by a factor of two when the signal-to-noise ratio is enhanced by a factor of $2$⁠. This was verified by injecting additive white Gaussian noise on top of the swept sine wave and noting that the threshold for lock loss at low τ reproduced this predicted behavior.

In terms of noise in the error signal, the most significant source is the product of noise in the input data since this noise enters the error signal without bandpassing through the IWAVE filter. However, this broadband noise has an rms amplitude independent of the amplitude of the tracked wave. If the tracked wave drops in amplitude but the rms of the broadband noise at the input does not, then the noise component of the error signal will scale inversely proportional to the amplitude of the line. This is exactly what is seen when IWAVE is run, in practice, on real data. The normalization of the error signal with the amplitude squared is necessary to ensure that the feedback loop gain is independent of the wave amplitude.

The other function of the error signal is to provide an indication of whether or not the PLL is locked. With a non-stationary rms, this is difficult. However, if we scale the error signal by multiplying it by the amplitude of the wave being tracked, this stabilizes the rms of the error signal noise at the level of the line amplitude. If we further divide by the long-term rms of the input data, we then obtain a statistic that has a stable rms of order unity when the servo is locked. Departures from lock manifest themselves as large transient spikes of amplitude greater than ten. This effect will be seen in the discussion of performance on gravitational wave data in Sec. IV. In particular, in Fig. 11, the error signal plotted for each of the four harmonics in the study has been scaled in this way.

Multiple IWAVE filters can be applied to multiple harmonic waves in a single data stream. However, when those harmonics are frequencies spaced closely together, this causes crosstalk between the different filters. All harmonics present in the data will enter every IWAVE instance through the error signal. To mitigate this effect, we employ a cross-subtraction scheme, as illustrated in Fig. 10.

The schematic shows only two IWAVE instances, but the technique generalizes to any number of filters. Assume that the two IWAVE filters are locked on line frequencies f₁ and f₂. The D and Q phase outputs from each IWAVE instance are fed into a phase shifter, which estimates the wave signal one sample in the future at the frequency of the wave tracked by the filter. This wave sample is stored until the next input data sample when it is subtracted from the input data to the other filter. In this way, the input data to each IWAVE filter is purged of the harmonic being followed by the other filter. This technique has been determined to successfully lock multiplets of up to 20 filters, leaving the frequency estimates from each filter free of oscillations at the difference between frequencies in the multiplet, the effect seen in the absence of the cross-subtraction technique. The technique works with IWAVE because the two quadrature outputs can be used to generate an output shifted through an arbitrary phase shift and because the outputs are at the same amplitude as the input wave onto which IWAVE is locked.

We present an example of the application of IWAVE to a set of harmonic noise components of gravitational wave data taken from the LIGO open data center web site.⁴⁰ The data were acquired by the LIGO Hanford interferometer on November 30th, 2016 and consists of 800 s of calibrated strain data from the Hanford interferometer during an 800 s lock stretch. The data were preprocessed with a fourth order Butterworth highpass filter at 30 Hz, followed by four third order Chebyshev type 2 bandpass filters between 5 and 300 Hz, applied in series. Finally, a fine adjustment to the spectrum was made using a single real pole at 10 Hz. The resulting data are dominated by the 20–80 Hz band and are approximately white in that range. It is not a requirement that the input data to IWAVE be whitened, but if a spectral feature is to be successfully tracked by IWAVE, its peak should rise above the noise floor in the surrounding background; whitening ensures that this is the case. Eight harmonic features were identified from a broad power spectrum and were tracked using eight parallel IWAVE instances. These originate from various instrumental sources present in the Hanford instrument at the time when the data were acquired. Violin mode harmonics have been studied in detail by Cumming et al.³¹ 

The results at four of the identified frequencies are shown in Fig. 11. For each harmonic, the frequency, amplitude, and scaled error signal (as described in Sec. III B) are displayed. Amplitudes are in dimensionless strain units, so one represents the 4 km length of the LIGO detector arms. Although the lines are in some cases almost degenerate in frequency, there is no evidence of beats between the reconstructed frequencies. The error signal is roughly static with approximately unit rms although the distribution is non-Gaussian because of the modulation of the input noise by the sinusoidal IWAVE output. A non-statistical transient fluctuation in the modified error signal in the 36.7 Hz line at around 360 s corresponds to a jump in the IWAVE reconstructed frequency, indicating that IWAVE momentarily lost lock when either the frequency of this sinusoid shifted rapidly or IWAVE jumped between two almost-degenerate harmonics. A second loss of lock can be seen in the 37.3 Hz data at around 30 s accompanying a sudden drop in the amplitude of the harmonic. The eight IWAVE instances all successfully tracked their target harmonics with the reconstructed error signal providing a useful performance indicator.

We have described IWAVE, a novel orthogonal state generator, resonant filter, and phase locked loop for the dynamic tracking of harmonic waves. The method has a single input parameter, its response time. The algorithm has a low computational load so that many harmonics can be tracked in real time using a single central processing unit (CPU) core. The ability to track multiple closely spaced harmonics means that the method lends itself well to applications where there are dense “forests” of harmonics, such as in LIGO violin mode clusters and communication applications. IWAVE has been applied to LIGO strain data and used to study the character of violin modes in ultralow loss fused silica suspensions. There are many possible applications of the IWAVE method. It is complementary to existing PLL algorithms that we have described in the  Appendix. We have supplied software implementations of IWAVE in C with MATLAB and Python wrappers to encourage the community to find other applications and uses.³⁸ 

The authors can see several directions in which IWAVE could be improved. The error signal is susceptible to broadband noise contamination, and a narrower band alternative would be of benefit in applications with very weak signals although narrowbanding will reduce the responsiveness of the method to frequency changes. There are also applications where feedback is not important, for example, the use of IWAVE for novel resonators that is promising for resonant detectors of weak signals in physics, such as those of dark matter axions.⁴¹ We look forward to seeing what the community finds to do with our harmonic tracking algorithm.

The members of the Sheffield gravitational wave research group acknowledge the support of the Science and Technology Facilities Council under Grant Nos. ST/V005693/1, ST/V001752/1, ST/V001744/1, ST/V001019/1, and ST/R000336/1. I.J.H. was supported by the Hollows Scientific Foundation. L.S. acknowledges the support of the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), Project No. CE170100004, the United States National Science Foundation, and the LIGO Laboratory. M.F. acknowledges the support of the Fonds de la Recherche Scientifique-FNRS, Belgium, under Grant No. 4.4501.19.

This research has made use of data, software, and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org/), a service of the LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. The LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the United States National Science Foundation and operates under cooperative Agreement No. PHY-1764464. Advanced LIGO was built under Grant No. PHY-0823459. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded through the European Gravitational Observatory (EGO) by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN), and the Dutch Nikhef with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, and Spain.

The authors have no conflicts to disclose.

The data that support the findings of this study are available at Ref. 38. The data for the performance tests discussed in Sec. IV can be obtained from the LIGO open data center web site.⁴⁰ 

The now often-used Generalized Integrator-Based PLL, SOGI-PLL, was introduced in 2006.²³ SOGI is an orthogonal state generator, producing in-phase and quadrature-phase copies of the input wave analogous to the IWAVE D and Q outputs. These outputs are mixed with two quadratures from a reference oscillator, leading to an error signal in the quadrature-phase output that indicates frequency differences between the SOGI output and the reference. The error signal is fed to a proportional/integral filter. The filter output is added to a frequency offset input ω_n, and the result is used to adjust the coefficients of the SOGI filter to reflect the frequency change as well as to control the reference oscillator. Essentially, the SOGI algorithm is coupled to a conventional phase locked loop with a reference oscillator.

The SOGI algorithm was developed for the field of power grid monitoring, where frequency changes are a small fraction of a nominal constant value, commonly 50, 60 Hz, or harmonics of these. It is not designed for large departures from the frequency set at the ω_n input. The SOGI orthogonal state generator is designed using two s-plane integration stages, shown in the tan SOGI block. These s-plane filters are transformed to the digital domain by, for example, using a Tustin algorithm. The SOGI method does not track the wave amplitude although this can be done in a separate circuit, assuming the frequency of the wave is known, using homodyne detection, for example. The SOGI orthogonal state generator does not have the same transfer functions as the IWAVE one, and this is not surprising given the different methods by which they are obtained although the denominators of the two transfer functions are identical as they both represent a resonant response to a harmonic drive.

The EPLL although apparently seeming to differ from quadrature signal generation-based PLLs, such as SOGI-PLL, is closely related thereto.²⁸ A conventional PLL is embedded within a second feedback servo that uses the integral of the quadrature phase of the PLL output, multiplied by the unintegrated quadrature phase to synthesize a 2ω signal, which can then be subtracted from the input data, compensating for the 2ω component of the phase detector output inside the PLL. Furthermore, the common amplitudes are also equal to half the amplitude of the sinusoid in the input data. In the EPLL structure, the input data are normalized to the PLL by this amplitude, thereby removing the amplitude dependence of the phase detector.⁴² 

Unlike IWAVE, there is a necessity for two servos to be locked at once and several numerical parameters that must be adjusted, including parameters necessary to implement the s-plane design digitally. EPLL is susceptible to interference from other waves present in the input and to DC offsets. The latter two issues are addressed by prefiltering the input to the EPLL. EPLL has been widely implemented because of its comparative simplicity, ability to track wave amplitude, and the availability of the code implementing the method.⁴³